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CHAPTER  I 


INTRODUCTION 

This  report  deals  with  the  digital  redesign  of  the  pitch  control  system 
of  an  unstable  semi-active  terminal  homing  missile  system.  It  also  concerns 
with  the  design  of  a  coupled  high-order  multivariable  control  system  and 
the  realization  of  a  multi-port  controller. 

In  Chapter  II  we  extend  the  newly  developed  dominant-data  matching 
method,  which  was  the  research  result  supported  by  the  U.  S.  Army  Missile 
R&D  Command  under  DAAK  40-78-C-0017 ,  to  design  the  digital  controller  for 
the  pitch  control  system  of  an  unstable  semi-active  terminal  homing  missile 
system.  In  Chapter  III,  we  introduce  a  new  and  simple  method  for  coupled 
high-order  multivariable  control  system  designs.  In  Chapter  IV,  we  develop 
a  new  block  canonical  form  of  a  matrix  transfer  function  for  possible 
realization  of  a  multi-port  controller  without  using  integrators. 

Other  new  findings  of  this  research  are  reported  in  the  appendix. 
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A  DOMINANT-DATA  MATCHING  METHOD  FOR  DIGITAL  FILTERS 
AND  DIGITAL  CONTROL  SYSTEMS  MODELING  AND  DESIGN 

L.  S.  Shieh1,  Y.  F.  Chang1,  and  R.  E.  Yates2 

ABSTRACT 

A  dominant-data  matching  method  is  presented  for  obtaining 
a  reduced-order  discrete-data  pulse-transfer  function  from  either 
a  high-order  continuous-data  transfer  function  or  a  high-order 
discrete-data  pulse-transf er  function,  and  for  identifying  the 
pulse-transfer  function  of  a  system  from  available  experimental 
time  and  frequency  response  data.  The  method  may  also  be  applied 
to  the  digital  control  systems  design  problem  with  various  sampling 
periods.  The  same  method  can  be  used  for  digital  filter  designs 
if  the  filter  specifications  obtained  by  this  method  are  viewed 
as  control  specifications.  The  discrete-data  system  has  the 
exact  dominant  characteristic  performance  of  the  original  con¬ 
tinuous-data  or  digital  system.  The  relaxation  of  the  sampling 
period  requirement  and  the  flexibility  of  our  new  method  facilitate 
the  practical  industrial  implementation  and  application. 


1L.  S'.  Shieh  and  Y.  F.  Chang  are  with  the  Department  of 
Electrical  Engineering,  University  of  Houston,  Houston,  Texas 
77004. 

2 

R.  E.  Yates  is  with  the  Guidance  and  Control  Directorate, 
•  U.S.  Army  Missile  Command,  Redstone  Arsenal,  Alabama  35809. 
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I.  INTRODUCTION 


Most  practical  industrial  circuits  and  control  systems  are 
continuous-time  systems  for  which  analog  filters  and  controllers 
are  employed  to  improve  performance.  The  recent  availability 
of  high  performance,  low  cost  microprocessors  and  associated 
electronics  has  led  to  replacement  of  many  continuous  systems 
with  systems  employing  digital  filters  and  controllers.  Many 
techniques  have  been  developed  for  digital  control  systems  design 
[l] -  [4] .  Among  them  the  u-domain  bilinear  transformation  is 
often  applied  to  design  industrial  digital  controllers.  However, 
this  method  is  graphical  and  involves  "cut-and-try"  procedures. 
Recently,  Kuo  [5]  and  others  developed  an  optimal  discrete-time 
data  matching  method  for  the  redesign  of  a  continuous-data  system. 
Constant  controllers  instead  of  dynamic  digital  controllers  are 
mainly  employed  in  these  designs.  As  a  result,  good  performances 
of  redesigned  systems  can  be  achieved  if  the  frequency  of  the 
input  signal  is  sufficiently  lower  than  the  sampling  frequency. 

As  an  alternate  to  Kuo's  time-domain  approach,  Rattan  and  Yeh 
[6]  have  given  an  elegant  frequency-domain  method  for  the  redesign 
of  continuous-data  systems.  The  method  of  weighted  least-squares 
complex-curve  fitting  due  to  Levy  [7]  and  Sanathanan  and  Koerner 
[8]  has  been  successfully  extended  in  the  z-domain  to  determine 
a  dynamic  digital  controller.  As  a  result  of  these  efforts, 
better  performance  of  redesigned  systems  can  be  achieved.  However, 
this  method  is  restricted  to  systems  whose  controllers  are  se¬ 
lected  in  such  a  way  that  the  linear  solution  of  the  unknown 
constants  in  the  controller  is  possible.  On  the  other  hand, 
if  both  feed-forward  and  feedback  dynamic  digital  controllers 
are  employed  in  the  design,  the  closed-loop  pulse- transfer  function 
may  have  nonlinear  coefficients  of  the  unknown  constants  of  the 
controllers.  Therefore,  the  linear  solution  to  their  method 
may  not  be  valid.  Furthermore,  most  often  design  goals  are  as¬ 
signed  by  using  a  mixture  of  time-domain,  frequency-domain  and 
complex-domain  control  specifications  [9]  rather  than  a  set  of 

frequency-response  requirements.  Therefore,  complex-curve  fitting 
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methods  may  not  be  applicable.  In  this  paper,  a  computer-aided 
method  is  proposed  for  matching  the  dominant  data  of  a  high-order 
continuous-data  system,  or  a  discrete-data  system,  with  the  dominant 
response  of  a  low-order  digital  system  replacement.  Also,  methods 
are  given  for  system  identification  and  digital  controller  design 
of  these  systems. 

II.  DOMINANT  DATA  AND  DOMINANT-DATA  MATCHING  METHOD 


The  characteristics  of  a  control  system  or  a  filter  are 
often  expressed  by  either  a  time-response  curve,  or  a  frequency- 
response  curve  or  a  set  of  poles  and  zeros  in  the  complex  plane, 
or  both.  The  quantitative  description  of  the  steady-state  behavior 
is  characterized  by  its  final  value  as  t-*-®  and  by  the  value 
of  the  steady-state  frequency  response  as  oj-*-0.  On  the  other 
hand,  the  quantitative  description  of  transient  behavior  is  rep¬ 
resented  by  its  time-domain  control  specification  [9]  (for  example, 
the  percentage  overshoot  and  the  rise  time)  and  by  the  frequency- 
domain  control  specifications  (for  example,  the  maximum  value 
of  the  closed-loop  frequency  response  and  the  bandwidth) .  These 
specifications  which  are  defined  for  control  systems  can  be  con¬ 
sidered  specifications  of  analog  or  digital  filters.  This  is 
because  a  digital  system  (or  a  discrete-data  system)  can  be  viewed 
as  a  continuous-data  system  in  the  frequency  domain  when  z  =  e-503*1, 
where  T  is  a  sampling  period. 

Some  empirical  observations  or  rules  of  thumb  due  to  Axelby 
[10]  that  link  the  specifications  of  the  continuous-data  systems 
in  both  the  time  and  frequency  domains  are  as  follows: 


1. 


**t  Mp  -  sin  4> 


(la) 


m 


M^:  Maximum  value  of  unit-step  response 

Mpj  Maximum  value  of  the  closed-loop  frequency  response 
m 


Phase  margin. 


2.  M  a  (lb) 

e  u) 

c 

M  :  Maximum  value  of  the  error  of  the  unit-ramp  function. 


wc:  Gain-crossover  frequency 


3.  “p  S  uc  (lc) 

Up!  The  peak  value  frequency  or  the  frequency  when  occurs. 

4.  Mfc  S  wc  (Id) 


Mfc:  Maximum  value  of  the  unit-impulse  response. 


5-  t  £ 

p  (jl) 

C  C 


t  :  The  peak  value  time  or  the  time  when  M.  occurs. 
P  t 


i.  «  1-8 

6*  fcv  "  ~ZT 
c 


t  :  The  time  when  the  maximum  error  of  the  ramp  function 
with  respect  to  its  input  occurs. 

T-  *c  3  !T  .  (19> 

c 

• 

t  :  The  time  when  M.  occurs, 
c  t 

Other  rules  of  thumb  according  to  Truxal  [11]  are: 

8.  t  wb  5  0.6"  to  0.9"  (lh) 

tf:  The  rise  time  or  the  time  required  for  the  response 

to  go  from  10  to  90  percent  of  its  final  value. 

a).  :  The  bandwidth  in  rad/sec. 


The  delay  time  or  the  time  required  to  reach  50 
percent  of  its  final  value. 

The  velocity  error  constant. 
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The  above  rules  can  be  verified  by  using  the  standard  second 
order  transfer  function: 


Y  (s) 
R(s) 


U>. 


n 


2  x  ,r  ,  2 

s  +  2^tons  +  U)n 


(2) 


where  R(s)  and  Y (s)  are  the  input  and  output  functions,  respec¬ 
tively,  £  is  the  damping  ratio,  and  is  the  undamped  natural 
angular  frequency.  The  zero  of  the  system  is  located  at  infinity, 
and  the  finite  poles  are  in  the  complex  plane.  The  £  and  wn 
as  used  herein  are  defined  as  the  complex-domain  control  specifi¬ 
cations  . 


Recently  Shieh  et  al.  [12]  have  studied  the  relationships 
between  the  complex-domain  specifications  (£  and  ton)  and  the 
time-domain  and  frequency-domain  specifications  by  using  a  more 
sophisticated  model 


Y(s) 

R(s) 


T  tO  S  U) 

n  +  n 


s2  +  2£u)  s  +  u>  2 
n  n 


TS*  +  1 

s*2  +  2£s*  +  1 
(3) 


where  s*  =  s/u>n  is  a  normalized  complex  variable  and  t  indicates 
the  location  of  a  finite  zero.  The  relationships  are  as  follows: 


(4a) 


Mfc  =  1  +  e  *  n  p(x2  -  2t£  +  1 
2.  Mp,  wp,  K,  o>n,  and  x 

“p  -  1  -  252  ifT.O 

Mp  =  1/  -  K2)  if  T  =  0 


(4b) 


(4c) 

(4d) 
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7. 


“b'  ^ '  wn'  and  T 


U), 


=  “n[( 


1  +  T  - 


0  +  4 


1  +  T  - 


2£2)2  +  lj 


-11/2 


(4n) 


The  above  analytical  expressions  can  be  plotted  by  letting 
con  =  1  and  £  the  variable.  Once  a  specification  is  defined, 
other  corresponding  specifications  can  then  be  determined.  In 
other  words,  if  time-domain  specifications  are  assigned,  the 
corresponding  frequency-domain  and  complex-domain  specifications 
can  be  determined.  The  frequency-response  data  or  the  equivalent 
data  obtained  from  Eqs.  (1)  and  (4)  at  u  ,  oj  ,  w  ,  w  ,  w.  ,  etc., 

XI  u  71  P  u 

are  the  dominant  data  as  used  in  this  paper.  Other  important 
data  will  be  the  frequency-response  at  s  =  ju)=jO  because  these 
characterize  steady-state  behavior. 


Our  new  dominant-data  matching  method  matches  the  above 
dominant  data  of  a  continuous-data  or  a  discrete-data  system 
to  those  of  the  newly  designed  or  modeled  discrete-data  system. 
The  steps  involved  are  as  follows: 


Step  1: 


Step  2: 


Step  3: 


Determine  a  set  of  dominant  frequency-response 
data  from  the  assigned  time-domain  and  complex- 
domain  specifications  by  using  the  rules  and  results 
in  Eqs.  (1)  and  (4) . 

Assume  a  fixed  configuration  digital  system  and 
controllers  with  unknown  constants.  Determine 
the  open-loop  and  the  overall  pulse-transf er  function 
of  the  system. 

Formulate  a  set  of  linear/nonlinear  equations 
by  matching  the  unknown  constants  of  the  pulse- 
transfer  function  and  the  assigned  dominant  data. 
Solve  the  equations  by  using  the  multidimensional 
Newton-Raphson  method  [13] ,  available  as  a  library 
computer  program  package  (called  the  Z  system) 
in  many  digital  computers  [14]  . 
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Step  4:  Estimate  initial  value  for  the  numerical  solution 

of  the  Newton-Raphson  method  by  constructing  a 
crude  pulse-transfer  function  which  can  be  obtained 
by  a  complex-curve  fitting  method. 

Step  5:  Compare  the  results  with  the  assigned  specifications. 

When  the  dominant  data  are  obtained  from  either  a  high-order 
continuous-data  transfer  function  or  a  high-order  discrete-data 
pulse-transf er  function  and  a  low-order  pulse- transfer  function 
is  required;  this  is  a  model  reduction  problem.  If  the  dominant 
data  are  determined  from  an  experimental  set  of  time  and  fre¬ 
quency  data  and  the  corresponding  pulse- transfer  function  is 
required,  this  is  the  identification  problem.  The  order  of  the 
identified  pulse- transfer  function  depends  on  the  number  of  dominant 
data  parameters  used.  Therefore,  the  identified  pulse-transfer 
function  could  be  the  reduced-order  model  of  the  original  high- 
order  system.  The  above  two  problems  can  be  considered  the  modeling 
problem.  When  the  design  goals  are  specified  by  a  set  of  dominant 
data  and  the  digital  controllers  with  unknown  constants  are  designed 
to  match  the  desired  dominant  data,  this  is  the  design  problem 
for  digital  control  systems.  Applications  of  the  above  new  method 
will  be  described  in  the  following  sections. 

III.  MODELING  A  REDUCED-ORDER  PULSE-TRANSFER  FUNCTION 


We  use  a  real  stabilized  pitch  control  system  of  a  semiactive 
terminal  homing  missile  [15]  as  an  illustrative  model  to  show 
that  the  characteristics  of  the  transient-state  response  of  a 
system  can  be  estimated  from  the  dominant  frequency-response 
data  and  the  applications  of  the  proposed  method  to  the  identifi¬ 
cation  and  model  reduction  problems.  A  block  diagram  of  the 
missile  system  is  shown  in  Fig.  1.  The  closed-loop  high-order 
transfer  function  is 


Y(s) 

R(s) 


Gc(s)Go(s) 

=  T  +  Gc(s)G0(s)Hg(s) 


Ge(s) 

1  +  Ge(s) 


<s) 


8 


(5a) 


where 


Gc(s) 


=  the  stabilization  filter 


1#  6  (25  +  X)  (l25 

+  1 1 
/ 

l 

[(l5o) 

I2  +  1 

(fes)-  *  ij  m 

12  +1 

f°*8l 

1200J 

Is  +  1] 

(5b) 


Gq(s) 


The  transfer  function  of  the  actuator  and  aerodynamics 
of  the  missile  system 


324332. 316  (s  +  0.1933)  (s  +  65)  (s  +  1500) 

=  s(s  -  2.921)  (s'+"3.175)  (s  +  87.9  +  j95.5)  (s  +  112.5)  (s  i 

(5c) 

Hg(s)  =  The  transfer  function  of  the  gyro  =  1  (5d) 

Ge(s)  =  Gc(s)Go(s)  =  The  unstable  open-loop  transfer  function 


of  the  existing  stabilized  system. 


(5e) 


The  closed-loop  transfer  function  Te(s)  becomes 


Te(s)  = 


in  q 

b^s  b-^s  4*  •  •  •  4*  s  4*  b-^Q 


a0si;L  +  a^s10  + 


+  a1Qs  +  a 
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(6) 


where 


a0  "  1 

a^  =  1.923554000  x  10' 
a2  =  9.316239040  x  105 
a3  =  2.976950696  x  10f 

a.  =  6.231675318  x  10 

4 

a5  =  9.360329977  x  10 
a6  =  9.749923212  x  10 
a?  =  6.667397031  x  10 
a8  =  2.420405431  x  10 


10 


12 


14 


16 


18 


a9  =  2.911920560  x  10 


a1Q  =  2.419047424  x  10 
an  =  8.802158509  x  10 


18 


19 


18 


b0  =  0 
bL  =  0 
b2  =  0 
b3  =  0 
b4  =  0 


5 

b6 

b7 

b8 

b9 

*10 


1.494523312 

X 

1011 

2.563396371 

X 

1014 

5.017212044 

X 

1016 

2.926344345 

X 

1018 

4.610004670 

X 

1019 

8.802158509 

X 

1018 

na 
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The  Nyquist  plots  of  Ge(s)  and  Gq(s)  are  shown  in  Fig.  2.  The 
dominant  data  of  Gg  (s)  are: 

1.  Real  and  imaginary  parts  of  Gg(s)  at  s  =  jw  =  jO  are 


Re  [Ge  { j 0 )]  =  -2.103817 

Im[Ge(jO)]  =  ® 


2.  Gain  margin  Gem  of  this  system  GgCjto^)  is 


Q  =  |  _ i _ 

em  Ge(j‘JJii) 


Re  [Ggtjw^)] 


1  i 

-1.51 


where  the  phase-crossover  frequency  to^  is  given  by 

a)  =  1.9  rad/sec,  such  that  /G  ( ju>  )  =  -180° 
n  /  e  tt 

The  equivalent  real  and  imaginary  parts  of  G^ju)^)  at  *  1.9 

rad/sec  are 


Re  fG  ( jw  )1  =  -1.507944 

e  it  J 

Im  [g  ( jo>  )1  =  -0.006490205 

L  e  J  tt  J 


3.  Phase  margin  4>em  of  the  system  Ge  ( ju> ) .  is 

<t»em  =  180°  +  /Ge(iMc)  5  5.7787°  (7f) 

where  the  gain  crossover  frequency  to  is  given  by  to  =  3.2  rad/sec 

c  c 

so  that, 

|Ge(jwc)l  =  1  (7g) 

Equivalent  real  and  imaginary  parts  of  Ge  ( jwc)  are 

Re[Ge(jioc)]  =  -0.5939143  (7h) 

Im[Ge(jtoc)]  =  -0.09547478  (7  i) 

It  is  required  to  determine  a  reduced-order  pulse-transf er  function 
such  that  the  characteristics  of  the  identified  discrete-data 
model  agree  as  closely  as  possible  with  those  of  the  high-order 
continuous-data  system. 


Let  the  required  overall  pulse- transfer  function  be 


Tr(z)  = 


where  the  open-loop  pulse-transfer  function  Gr(z)  is 


.  +  x_  ,  z  +  x_ 
_ P"1  P 


XnZ^  +  X, zp  ^  +  . 

Gr(z)  =  - y - —i - — 3- 

(z  -  1) (yQzq  +  yxz4  +  ...  +  y^z  +  yq) 

Gr(z)  is  assigned  to  be  a  type  "1"  system  because  Gg  (s)  in  Eq. 
(5e)  is  a  type  "1"  system.  To  match  the  five  dominant  data  in 
Eq.  (7)  we  choose  q  =  2,  p  =  2,  and  y^  =  1.  Thus  Gr  (z)  becomes 


(8b) 


G  fz)  = 


XQZ  +  X  x  Z  +  x2 

(z  -  1)  (z2  +  yxz  +  y2) 


(9) 


where  x^  and  y^  in  Eq.  (9)  are  unknown  constants  to  be  determined. 
The  goal  is  to  determine  the  unknown  constants  x  and  y  in  G  (z) 
so  that  G  (z)  as  z  =  eJ  matches  the  dominant  data  in  Eq.  (7) . 

The  sampling  period  T  (=0.008  sec)  and  to_  (  =  250it  rad/sec)  are 
chosen  to  be  synchronized  with  the  125  Hz  pulse-width  modulated 
actuator  [16].  Since  Ge(s)  and  Gr(z)  are  type  "1"  systems  and 
we  need  to  match  the  dominant  data  of  Ge(jw)  as  <0  =  0  in  Eq. 

(7a)  and  the  dominant  data  of  Gr(z)  as  (a)  =  0  or  z  =  e-*uT  =  1 
in  Eqs.  (8b)  or  (9),  the  expression  for  Gr(z)  in  Eq.  (8b)  is 
modified  as  follows: 


Substituting  z  =  z*  +  1  into  Eq.  (8b)  yields 


Gr(z*)  = 


*  *  *  *  *n 

XP  +  Viz  +  +  xoz 


*  *  *  ★  *Q 

2*(yq  +  yq-lZ  +  +  y0Z  ] 


*-l  *  *2 

=  e_1z  +  eQ  +  e^z  +  e2z  +  ... 


(10a) 


where 


p  p  q 

xp  =  J0xi'  xp-i  =  ifi^p-i'  yq  =  i=0yi'  yq-i 


=  ..E.iyq-i' 


i=l 


e-l  =  Xp/yq  nnd  C0  =  (yqXp-l  “  yq-lXp)/yq2'  etC* 


Equating  the  respective  real  and  imaginary  parts  of  Gr(z)  for 
to  =  0  and  those  of  Gr  (z*)  for  to  =  0  gives 
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K(z)l  =  Re  [g  (z*  j] 

L  r  J  z  =  1  +  jO  L  r  J  z*  =  z-1  =  jO 

[Gr(z)l  =  lm[Gr(z*)l 

L  1  J  z  =  1  +  jO  L  E  J  z*  =  z-1  =  iO 


(10b) 


(10c) 


Eqs.  (10b)  and  (10c)  imply  that  eQ  in  Eq.  (10b)  is  the  asymptotic 
line  of  the  type  "1"  systems  at  low  frequencies. 

In  the  frequency  domain,  Eq.  (9)  can  be  expressed  in  an 
alternative  form  as  follows: 

Let  us  define 


3U|rT  A 

z  =  e  =  cos  <i)kT  +  j  sin  wkT  =  uk  +  jvk 
and  substituting  z  =  uk  +  jvk  into  Eq.  (9) ,  we  have 


(Ha) 


2  2 

(XqUj^  ^  x^uk  +  X2)  j  (2XgUkvk  +  x^Vj^) 

k'  k  [,uk  -1)  (uk2  -  vk  +  yiuk  +  y2)  -  vk(2ukvk  +  y^)] 


+  j  [(uk  -l)  (2ukvk  +  yxvk)  +  vk(uk2  -  vk2  +  Yluk  +  y2)] 


=  R.  +  j  I, 


(lib) 


where  wk  are  specific  frequencies  and  Rk  =  Re  jGr  (uk , vk )J  ,  Ik  = 

Im  |Gr  (uk , vk)J  .  If  Rk  and  Ik  are  the  known  or  assigned  values  at 
frequencies  wk,  we  can  obtain  two  linear  equations.  First,  we 
multiply  both  sides  of  Eq.  (11b)  by  the  common  denominator,  then 
we  separate  the  real  and  imaginary  parts  and  then  equate  the 
respective  real  and  imaginary  parts.  Thus  we  have 

2  2 

f i  (xo,xl'x2'yl'y2*  =  (x0uk  “  xovk  +  xluk  +  X2J 

-  Rk  [(uk  -  1)  (uk2  -  vk2  +  ylU|i  +  y2) 

-  V2ukvk  +  ^lV]  +  ^['“k  '  1><2ukvk 

+  y^)  +  vk(uk2  -  vk2  +  +  y2)]  =  0 

12  die) 


and 


fi+l(x0'xl'x2'yl,y2)  =  <2x0ukvk  +  xlvk) 

-  Rk  [(uk  - X)  (2ukvk +  yiV  +  vuk2  -  vk2 
+  yiuk +  y2>]  -  h  K  - 1}  {uk2  -  vk2 
+  +  y2)  -  vk(2ukvk  +  yivk}]  =  0 

(lid) 

Using  the  expressions  in  Eqs.  (10b),  (11c),  and  (lid)  and  the 
assigned  dominant  data  in  Eq.  (7)  we  can  formulate  one  nonlinear 
equation  and  four  linear  equations  ^(Xj/y^)  =  0  for  1  =  1/2,...,  5 
as  follows: 

(i)  The  data  in  Eq.  (7a),  or  Re  [Ge  ( jm)J  =  -2.103  for  0 
and  the  relationship  in  Eq.  (10b)  gives  a  nonlinear  equation: 

f 1 <Xfi.  )  =  (2x0  +  Xl}  ^  +  yl  +  y2*  "  (2  +  (xo  +  xi  +  x2> 

+  2.103(1  +  yx  +  y2)2  =  0  (12) 

(ii)  From  Eqs.  (7d)  and  (7e)  ,  Re  [Ge  (3“^ )]  =  -1.507944  and 

Im  [Ge  ( ju^)]  =  -0.00649025  for  =  1.9  rad/sec.  We  define 

RR  =  Ri  9  =  -1.507944,  Ik  =  ^  g  =  -0.006490205,  u)R  =  “i  g  =  %  = 

uk  =  ui  9  =  cos  0)-^  gT  =  0.99988448  and  vk  ^  vi  9  =  sin  “1  9T 

=  0.01519941  as  T  =  0.008  sec 

Substituting  the  above  data  into  Eqs.  (11c)  and  (lid)  and  letting 
i  =  2  we  have  tv/o  linear  equations  f2(x^,y^)  =  0  and  ^(Xj^y^)  =  0 
as  shown  in  Eqs.  (11c)  and  (lid). 

(iii)  From  Eqs.  ( 7 h )  and  (7i)  we  define  u>k  =  coc  =  3.2  =  2» 

uk  =  u3  2  =  cos  w3  2T  =  °-99967234'  vk  "  v3.2  =  sin  w3  2T  =  °-025597204 
Rk  “  R3  2  =  -0 •  "39143  and  "  x3  2  “  -0*0"47478.  Substituting 

the  above  data  and  i  =  4  into  Eqs.  (11c)  and  (lid)  yields  two 

more  linear  equations  f4(x^,y^)  =0  and  fg(x^,ya)  =  0.  Thus, 
we  have  five  simultaneous  equations  f^x^y^)  =  0  with  five  unknown 
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constants  and  y^  to  be  solved.  Notice  that  if  the  data  of 
Eqs.  (7b) ,  (7c) ,  ( 7 f ) ,  and  (7g)  are  used  to  match  the  unknown 
coefficients  of  Eq.  (lib),  the  resulting  equations  f.(x0,y.)  =  0 
in  general  are  nonlinear.  Therefore,  we  note  that  in  general 
the  equations  f  ^  )  =  0  are  nonlinear.  The  Newton-Raphson 

method  available  as  a  library  computer  program  in  most  digital 
computers  [14]  can  be  applied  to  solve  these  nonlinear  equations. 
However,  as  is  well  known,  the  Newton-Raphson  method  will  converge 
to  a  desired  solution  for  a  small  range  of  starting  values  or 
initial  solution  estimates.  To  improve  the  convergence  and  to 
obtain  the  set  of  desired  solutions,  we  offer  the  following  method 
for  initial  estimates. 


Since  ^(x^y^)  =  0  is  nonlinear  and 


Yo)  “  0,  i  =  2 , .  . 


are  linear  equations,  we  linearize  ^(x^y^)  =  0  by  choosing 

a  very  low  frequency.  For  example,  if  we  choose  oj^  =  0.01  = 
^  n  _  *1  1  ▼  A—  .  /-i  ■*  — l  A 


W0.01'  then  Rk  "  R0.01 


0.01 


=  40.17319,  u.  =  u 


=  -2.1,  Ik  = 

COS  u0.0lV  0-99999950  and  vR  £  vQ>01  =  sin  a)0>0]T  =  8 


x  10 


o5oi 


Solving  ^  (x£,ya)  =  0  and  ^(x^y^)  =  0,  i  =  2,. ..,5  for  the 
unknown  constants  x^  (defined  as  x^*)  and  y&  (defined  as  y^*) 
we  get  xQ*  =  0.00679254,  x^  =  -0.0123359,  x2*  =  0.00554537, 
y^*  =  -1.9985417,  and  y2*  =  0.99794573.  Using  these  values  as 
initial  estimates  for  the  solutions  of  f  j_  (x^/Y^)  =  0  using  the 
Newton-Raphson  method  we  obtain  the  solution  x^  =0.00679259, 
xx  =  -0.01233599,  x2  =  0.00554531,  yL  =  -1.9985412,  and  y2  = 
0.9979452  at  the  second  iteration  with  error  tolerance  of  10-^. 
The  .desired  open-loop  pulse- transfer  function  is 


;r(z)  = 


0 . 006792596z  -  0.012335992z  +  0.0055453114 

z3  -  2 . 99854 12 z3  +  2.9964864z  -  0.9979452 


(13) 


A  Nyquist  plot  of  Gf(z)  is  shown  in  Fig.  2.  The  plot  matches 
closely  that  of  GR(s)  not  only  at  the  dominant  frequencies  but 
also  at  others.  The  Gr(z)  is  seen  to  be  a  good  reduced  model 
of  the  original  unstable  system  G  (s) .  This  is  the  contribution 
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of  our  new  method  because  there  are  no  known  effective  model- 
reduction  methods  for  unstable  systems.  The  resulting  closed- 
loop  pulse-transfer  function  which  is  the  reduced-order  discrete- 
data  model  of  the  original  high-order  continuous  data  system 
is 


Gr  (Z)  =  0.006792596Z2  -  0.0 12335992z  +  0.0055453114 

+  Gr(z)  z 3  _  2. 991748524z2  +  2.984150408Z  -  0.9923998886 


(14) 

Since  the  assigned  dominant  data  are  the  steady-state  frequency 
response,  it  is  interesting  to  compare  responses  of  Tg  (s)  in 
Eq.  (6)  and  Tf(z)  in  Eq.  (14)  shown  in  Fig.  3.  Observe  that 
both  the  transient  response  and  steady-state  response  of  the 
reduced-order  model  Tr(z)  are  excellent  matches  of  the  original 
high-order  system.  This  indicates  that  the  dynamic  characteristics 
of  the  system  (for  example,  peak  value  time  and  overshoot,  which 
may  not  occur  at  the  sampling  time)  are  indirectly  controlled 
by  the  assignment  of  the  gain-crossover  frequency  and  the  phase 
margin.  This  is  a  major  advantage  of  our  new  method.  Also  note 
that  the  reduced-order  model  gives  an  excellent  approximation 
of  the  original  system  when  driven  by  high-frequency  input  signals. 

To  determine  the  initial  estimates  x^*  and  y^*,  a  general 
formulation  of  a  set  of  linear  equations  can  be  constructed  from 
the  following  complex-curve  fitting  method. 


Consider  the  pulse- transfer  function 


*  m 
x0z 

+ 

x*zra_1 

+ 

...  + 

* 

X 

m 

*  n 

y0z 

+ 

+ 

...  + 

* 

*n 

(15a) 


*  *  * 

where  yn  =  1  and  x0  and  y0  are  unknown  constants  to  be  determined. 
^  r  ^  rii)  T  * 

Substituting  z  =  eJ  k  =  cos  rw^T  +  j  sin  rw^T  into  Eq.  (15a) 
gives 
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(15b) 


G 


111  ^  111  ^ 

Z  x  cos(m  -  Dui.T  +  j  I  x  sin  (m  -  S,)u.T 

_  g,  =  0  _  1  =  0  1 _ _ 

n  *  n  * 

Z  y.  cos  (n  -  +  j  I  y.sin(n  -  J,)tu.T 

S,=0  *  *  £  =  0  K 


=  R(tok)  +  jl(wk)  =  +  jl^ 

where  R^  and  1^  are  the  real  and  imaginary  parts  of  the  transfer 
function  at  the  experimental  frequencies  or  the  assigned  frequencies 
wk.  After  multiplying  both  sides  of  Eq.  (15b)  by  the  common 
denominator  and  separating  the  real  and  imaginary  parts,  we  equate 
the  respective  real  and  imaginary  parts.  This  yields  the  following 


matrix 

equation : 

COSKkj^T 

cos(m-lJw^T 

1 

(-R^cos (n-1 ♦ 

l1sin|n-J)u1T) 

(-R^cosuj^T  ♦ 

I^sinwjT) 

-R1 

sinmu^T 

sin{n-l)u1T 

0 

(-R^in  (n-lJw^T  - 

IjCOSln-llu^T) 

(-RjSinu^T  - 

^cosaijT) 

cosmwjT 

coslm-lluijT 

1 

(-R^cos (n-l)u^T  ♦ 

Ij*in (n-l)wjT) 

(-R^cosujT  ♦ 

IjSinujT) 

-Ri 

sinmujT 

ainlm-DujT 

0 

(-RAsin  (n-lju^  - 

I jCOs (n-l)ujT) 

(-RjSinu^T  - 

IicosuiT) 

t 

( 

i 


(RgCOsnuigT  -  I0Blnnw0T) 
(R0slnnu0T  ♦  IgCOsnugT) 


■  (RjCOanujT  -  IjSinnujT) 
(RjSlnnidjT  ♦  IjCosnuijT) 


(15c) 


Substituting  the  selected  (n  +  m  +  1)  frequency-response  data 

into  Eq.  (15c)  ,  we  can  solve  for  the  required  (n  +  m  +  1)  unknown 
*  .  * 

constants  x.  and  y  . 

^  **  *1  £ 


IV.  DIGITAL  CONTROL  SYSTEM  DESIGN 


Consider  the  pitch  control  transfer  function  of  the  missile 
system  of  Eq.  (5a).  The  unity-feedback  system  without  the  sta¬ 
bilization  filter  Gc(s)  is  unstable,  and  a  rate  gyro  is  not  available 
for  this  example  system.  It  is  required  to  design  a  digital 

controller  G  (z)  instead  of  an  analog  controller  G  (s)  such 
c  c 

that  the  designed  system  has  the  exact  control  specifications 
[9]  of  the  original  stabilized  continuous-data  system  given  in 
Eq.  (7).  Furthermore,  the  response  Gg  { jto)  at  w  =  140  rad/sec  = 

“l'lO  chosen  as  a  dominant  data  constraint  because  the  system 
has  an  inherent  high  frequency  noise  component  at  This 

is  a  digital  redesign  problem.  The  structure  of  the  digital 
control  system  is  shown  in  Fig.  4.  The  closed-loop  pulse- transfer 
function  of  the  desired  digital  system  becomes 


Y  (z)  _  Gc(z)GnGo(z) 


Gel<z> 


R(z) 


1  +  Gc(z)GnG0(z)  =  Tel(z)  =  1  +  G~[^) 


where 


G  G  (z)  =  (1  -  z~1)Z[-  G_ 
no''  |_s  o 

(.»]  ■  b°zv  +  bl\  +  •  •  •  + 

aQz  +  a^z  +  ...  +  a? 

bQ  =  0.4095517916  x  10”4 

a0  =  1'° 

bx  =  0.2526111734  x  10~3 

a2  =  -4.120000127 

b2  =  -0.2575534058  x  10”3 

a2  =  6.894911119 

b3  =  -0.1540502340  x  10-3 

a3  =  -6.064285805 

b4  =  0.1096039580  x  10"3 

a4  =  3.023078996 

b5  =  0.8808224830  x  10_5 

a5  =  -0.83313251 

b6  =  -0.1358278543  x  10-9 

a6  =  0.09942985948 

a?  =  -0.1532773068  x  10~5 

elU)  =  Gr(z)GnGJz) 
1  -  e 


c  '  -n~o 
-sT 


Gn(s)  = 


=  the  zero-order  hold 


(16) 


T  =  the  sampling  period  =  0.008  sec. 
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T 


G_(z)  is  the  desired  digital  controller  and  is 

G 


XnZ  +  X  2  4  X, 

Gc(z,  =  ± 2 

z  +  y.z  4-  y 


(17) 


where  x^  and  y^  are  unknown  constants  to  be  determined. 
G  (z)  is  a  forward  controller,  the  equations  f.(x.,y  ) 
be  formulated  from  the  following  equations: 


Ge(s) 


gc(2) 


s  =  gw. 


z  =  e3tokT 


VoWlz  =  e^kT 


Because 
0  can 


(18) 


where  o)q  =  0,  g  =  1.9,  “3  2  =  3.2,  and  =  140.  Notice 

that  G  G  (z)  for  z  =  e3tokT  is  not  equal  to  G  (s)  for  s  =  iio. 
no  o  J  k 

unless  T  =  0.  Using  the  dominant  data  of  Eq.  (7),  the  required 
response  at  cj140,  and  the  relationships  expressed  in  Eqs.  (10b) 
and  (18)  yields  a  set  of  equations  =  0  for  i  =  1,2,..., 5 

as  follows: 


(i)  Using  Eqs.  (7a)  and  (10b)  when  w  =  Wq  =  0,  yields  a 
nonlinear  equation; 

fi(xi,yJl)  =  -4.557577105  x  10~8 (xQ  +  x±  +  x2)  (1  4-  y;L  4-  y2) 

-  7.016133905  x  10-11  [(2xQ  4-  x±)  (1  +  y1  +  y2) 

“  (2  4-  yx)  (xQ  +  x1  +  x2)]  4-  2.1  x  3.505134712  x  10~8 

(1  +  y1  +  y2) 2  =  0  (19a) 

(ii)  Using  Eqs.  (7 h)  ,  (7i),  and  (18)  when  w  =  to3  2  =  3.2 
we  ge  t : 

r  2^*  1  r  3u3  2 ^  ~\ 

Re  Ig  (e  )J  =  1.5987861  and  Im|Gc(e  )J  =  0.23560917 


The  resulting  linear  equation  is 

2  2 

^2  ^Jl^Jl)  ~  ( xqu3  2  xqv3  2  xiu3  2  x2^ 


-  1.5987861(u3i22  -  v3_22  +  y3u3_2  +  y2) 


4-  0.23560917  (2u3  2v3.2  +  ylv3.2}  =  0 
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(19b) 


where  u3  =  0.99967234  and  v3  2  =  0.025597204 

(iii)  Using  Eqs.  (18)  and  (19b)  for  oj  =  w3  2  =  3.2,  obtain 
a  linear  equation 

f  3  ( » y )  =  (2xou3  2V3  2  +  xlv3  2*  ~  1* 5987861  (2u3  2V3  2 

+  Ylv3.2>  -  0-235609n(u3_22  -  v322 
+  y1u3<2  +  y2)  =  0  (19c) 

(iv)  From  Eq.  (18)  for  oj  =  uj^4q  =  140  we  have 

r  3U)140T  1  T  3t0l40Tl 

Re|Gc(e  )J  =  26.951878  and  Im  [ReGc  (e  )J  =  19.196865 

The  resulting  linear  equation  is 

f4(X«.,y)l)  =  (X0U140  "  X0V140  +  xlu140  +  X2) 

26.951878 (u14Q  -  v14Q  +  y^u^g  +  y2) 

+  19.196865 (2u140v140  +  y1v14Q)  =  0  (19d) 

where  u^4g  =  0.43568245  and  v^4g  =  0.90010044. 

(v)  From  Eq.  (19d)  for  w  =  u14Qf  we  obtain  another  linear 
equation 

f5(W  =  (2x0U140  +  Xlv140^  ~  26  -951878  (2ui40vi40  +  ylv140) 

-  19.196865  (u14Q2  -  v14Q2  +  y^g  +  y2)  =  0  (19e) 

The  above  set  of  linear  and  nonlinear  equations  can  be  solved 

using  the  Newton-Raphson  method.  The  initial  estimates  for  the 

Newton-Raphson  solution  may  be  determined  from  Eq.  (15c) .  Another 
,  * 

linear  equation  f^(x^,y^)  =  0,  instead  of  f  ^  ( x£ '  y^  ^  =  9  in  Eq. 

(19a) ,  can  be  constructed  to  yield  five  linear  equations  with 

*  *  A 

five  unknown  constants  (x^  and  y  ).  Ge(}a>)  at  u>  =  0.01  =  ojg 

is  used  in  this  case.  Substituting  Re  [g  (e7a)0 . 01T)]  =  1.5961120 

and  Im  [bc  (e7u0 . 01T)]  =  6.409642  into  Eq.  (18)  we  get  the  linear 

equation  19 


fl(Xi'yil)  =  (x0u0.01  ”  X0v0.01  +  Xlu0.01  +  X2} 

-  1.5961120  (u0>012  -  v0>01  +  +  y2) 

+  6.409642  (2u0  01v0  01  +  =  0  (20) 

where  uQ  =  0.99999950  and  vq  oi  =  8  x  10  5. 

* 

Solving  fi(xi/Yji)  *  0,  i  =  2,.. .,5  in  Eq.  (19)  and  f]_(x4*yA)  =  0 
in  Eq.  (20)  gives  a  set  of  initial  values  x^  and  y^;  x^  =  2.9025918, 
x*  =  9.358033b,  x*  =  -10.066413,  y*  =  -0.42702018,  and  y*  =  0.80224909. 
Using  these  initial  values  and  the  Newton-Raphson  method  to  solve 
the  equations  in  Eq.  (19)  we  obtain  the  solution:  Xq  =  11.869083, 

Xjl  =  -13.49237,  x2  =  3.0584008  ,  y1  -  -0.75055299  ,  and  y2  =  0.64699237 
at  the  fifth  iteration  with  the  error  tolerance  of  10-6.  The 
required  digital  compensator  is 

„  .  .  11 . 869 08 3 z2  -  13 . 49  237 z  +  3.0584008 

G  (Z)  =  - - - 

c  z*  -  0 . 75055299 z  +  0.64699237 

11. 869083  (z  -  0.82408067)  (z  -  0.31268533) 

=  (z  -  0.37527650  +  jO. 71144917) (z  -  0.37527650  -  jO. 71144917) 

(21a) 

The  Nyquist  plot  of  Gel(z)  =  Gc(z)GnG0(z)  is  shown  in  Fig.  2. 

It  closely  matches  the  Nyquist  plot  of  Gg  (s) .  The  closed-loop 
pulse- transfer  function  is 

T  =  0.4861004207xl0~3z8  +  0 ; 24 4 56 80553xl0~2 z7  -  0 . 63399888 15xl0~2z6 

el  z9  -  4 . 870067017  z8  +  10 . 63662758 z7  -  13.9112306z6  +  12.03802088z 

+  0.2419157047xl0~2z5  +  0 , 259 169 9 688 x!0~2 z4  -  0 . 1845418962xl0~2z 
-  7 . 023068 4  35  z^  +  2 . 678803581z3  -  0 . 6134 429209 z2 

+  0.2163673922xl0~3z2  +  0 . 269409 1451xl0~4 z  -  0 . 4 154 160 183xl0~9 
+  0.6435845178xl0-1z  -  0 . 9921078960xl0-6 
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(21b) 


The  unit-step  responses  of  the  existing  stabilized  continuous- 
data  system  Te(s)  in  Eq.  (6)  and  Tel(z)  in  Eq.  (21b)  are  shown 
for  comparison  in  Fig.  3.  The  time-response  of  the  newly  designed 
sample  data  system  is  very  close  to  the  existing  stabilized  system. 
It  is  interesting  to  note  that  Gc(s)  of  Eq.  (5b)  is  a  fourth- 
order  analog  controller  whereas  the  Gc(z)  of  Eq.  (21a)  is  a  second- 
order  digital  controller. 


In  a  large  control  system  it  is  often  difficult  to  select 
a  minimum  common  sampling  period  among  the  various  subsystems. 

For  example,  the  missile  inner  loop  stability  system  with  sampling 
period  T  =  0.008  sec  is  used  with  a  terminal  guidance  system. 

The  terminal  guidance  system  is  low-pass.  Therefore,  a  larger 
sampling  period  may  be  economically  used  in  this  system.  If 
we  assign  a  larger  sampling  period  T^  (  =  0.015  sec)  for  the  outer 
guidance  loop,  and  we  desire  a  single  sample  period,  we  must  raise 
the  sampling  period  T  (  =  0.008  sec)  of  the  actuator  and  inner  loop 
from  T  ( =0 . 008 )  to  T^  (  =  0.015).  Notice  that  the  new  sampling  fre¬ 
quency  cosl  (  =  2tt/T-^  =  418 . 88 )  >>  2a)140  (  =  280  rad/sec).  The  modified 
open-loop  pulse- transfer  function  with  Tj^  =  0.015  sec  is 


GnGo<z> 


C.  C 

SqZ  +  a^z  +  ...  +  ag 

—  —  —  - 

bQz  +  b1z  +  ...  +  b? 


where 

aQ  =  0.3733134407  x  10" 
a1  =  0.1570750710  x  10" 
a2  =  -0.169231262  x  10" 
a3  =  -0.5827024448  x  10 
a4  =  0.3184265948  x  10" 
a5  =  0.1895463329  x  10" 
a6  =  0.11354062  x  10-9 


(22) 


1.0 

-3.257024652 

3.855486034 

-2.039756267 

0.5526488259 

-0.1245437594 

0.1318981958  x  10“ 

-0.1252480548  x  10 


X 
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Since  a  larger  sampling  period  T.  is  used,  we  select  a  third 

* 

order  digital  controller  Gc(z)  rather  than  the  second  order  digital 
controller.  Gc(z)  of  Eq.  (21a)  is 


3  2 

*  xnz  +  x.z  +  x_z  +  x, 

G(z)  =  - h - - - "  (23) 

z  +  +  y2z  +  y3 

The  x0  and  y.  are  seven  unknown  constants  to  be  determined. 

^  A  A 

Ge(ju>)  at  a)  =  0  =  u)q,  u)  =  1.9  =  ui^  g,  co  =  3.2  =  u>3  2,  an^  **>  =  140 

u)i40  shown  in  Eq.  (7)  are  used  as  the  dominant  data  to  determine 

x^  and  y^.  Using  the  above  design  procedure,  we  can  determine 

a  set  of  equations  =  0,  i  =  1,2,.. .,7.  These  equations 

can  be  solved  by  using  the  Newton-Raphson  method,  with  the  set 

of  initial  estimates  obtained  from  Eq.  (15c) .  The  data  obtained 

from  Eq.  (18)  at  m  =  0.01  =  ojg  gX,  g,  2,  and  “^.40  are  use<^ 

in  Eq.  (15c)  to  determine  the  initial  estimates  Xg  =  13.177031, 

x*  =  -25.535836,  x*  =  14.643983,  x*  =  -2.2787512,  y*  =  -0.37332775, 

y*  =  -0.32614954,  and  y*  =  -0.296499004.  Using  these  values 
as  initial  values  for  the  Newton-Raphson  method,  gives  the  desired 
constants  x  and  y.  at  the  17th  iteration  with  error  tolerance 
of  10  .  The  newly  designed  digital  controller  G  (z)  is 

„*  ,  .  13 . 170704 z3  -  25.531430Z2  +  14.629635z  -  2.2685451 

G  (z)  =  - - - - - 

c  z  -  0 . 37424841z^  -  0.32757047z  -  0.29794756 

(24) 


The 


closed-loop  pulse-transf er 


V(z) 

R(z) 


=  T 


e2 


(z) 


Ge2(Z) 

1  +  Ge2(z) 


function  is 

bg  z9  +  b^  z8  + 
_  _ 

a  g  z  +  a^z  + 


+  b. 


+  a 


10 


(25) 


where 

Ge2(z)  =  G*  (  z)  G^G*  ( z) 

aQ  -  1.0 

ax  =  -3.626356261 
a2  =  4.758008528 


bg  =  0.4916800827  x  10 
bx  =  0.1115666668  x  10-1 
b2  =  -0.5693102109  x  10~ 
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<n 


a3  =  -2.77063927 

a.  =  1.081168732 

4 

a5  =  -0.8211905469 
a6  =  0.4739432136 
a?  =  -0.1233033662 
ag  =  0.3234184521  x  10" 
a9  =  -0.3972872336  x  10 
a1Q  =  -0.253840284  x  10~ 


b3  =  0.5766519111  x  10 
b4  =  -0.9250105745  x  10 
b5  =  -0.1256587702  x  10 
bg  =  0.5496414235  x  10~ 
b?  =  -0.4450686236  x  10 
bg  =  -0.4299777941  x  10 
bg  =  -0.2575720172  x  10 


The  Nyquist  plot  of  Ge2(z)  shown  in  Fig.  2  matches  very  well 
that  of  G  (s) .  The  unit-step  response  for  Y(z)  in  Eq.  (25)  is 
shown  in  Fig.  3.  The  time  response  of  Te2(z)  very  closely  matches 
that  of  the  original  system  Te(s).  The  resulting  design  is  seen 
to  be  quite  satisfactory. 


V.  CONCLUSION 

A  dominant-data  matching  method  has  been  given  for  fitting 
the  coefficients  of  a  pulse-transf er  function  from  available 
time  and  frequency  response  data  or  from  assigned  design  goals 
expressed  by  a  set  of  control  specifications.  When  the  dominant 
data  are  obtained  fron  a  high-order  continuous-data  as  well  as 
a  discrete-data  system,  our  new  method  has  been  used  to  determine 
a  reduced-order  discrete-data  system.  If  the  data  are  experi¬ 
mental  time  and  frequency  response  data  of  a  system  to  be  identi¬ 
fied,  our  method  may  be  used  to  identify  the  pulse- transfer  function. 
Also,  the  method  has  been  used  for  redesigns  of  a  continuous 
system  using  a  digital  filter  with  various  sampling  periods. 

The  pulse- transf er  function  obtained  by  our  new  method  has  the 
exact  dominant  performance  of  the  original  or  desired  system. 

We  feel  that  the  flexibility  and  accuracy  of  our  new  method  will 
have  significant  practical  advantages  for  the  design  of  digital 
systems . 
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igure  1.  Block  Diagram  of  a  Missile  Pitch  Control  System 


Figure  4.  Block  Diagram  of  a  Digital  Pitch  Control  System 
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Chapter  l ll 


A  MODIFIED  DIRECT -DECOUP LING  METHOD  FOR 
MULTIVARIABLE  CONTROL  SYSTEM  DESIGNS 

L.  S.  Shieh1,  Y.  J.  Wei1,  and  R.  E.  Yates2 

ABSTRACT 

A  design  method,  which  decouples  an  interactive  system  by 
using  a  compensator  obtained  from  the  plant  inverse  matrix,  which 
is  often  called  the  direct-decoupling  method  is  modified  in  this 
paper.  The  modified  direct-decoupling  method  uses  the  adjoint 
matrix  instead  of  the  inverse  of  the  plant  matrix  to  construct 
the  compensator.  The  method  uses  a  frequency-domain  model- 
reduction  method  to  simplify  the  degree  of  the  given  plant  transfer 
function  matrix  and  the  obtained  compensator.  For  an  open-loop 
stable  multivariable  system,  the  proposed  method  gives  a  simple, 
practical  and  realizable  controller  without  using  an  unstable 
pole-zero  cancellation  approach. 


1L.  S.  Shieh  and  Y.  J.  Wei  are  with  the  Department  of  Elec¬ 
trical  Engineering,  University  of  Houston,  Houston,  Texas  77004. 
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Arsenal,  Alabama  35809. 


1 


X.  INTRODUCTION 


For  the  general  multi-input-output  control  system,  each 
input  affects  several  outputs  and  there  are  many  degrees  of  free¬ 
dom  for  system  design.  Therefore,  it  is  difficult  to  control 
and/or  design  such  a  system.  The  removal  of  the  interactive  con¬ 
trol  effects  and  the  application  of  the  well-developed  classical 
design  techniques  of  a  single-variable  system  to  the  decoupled 
system  is  one  popular  design  method.  This  method  is  often  called 
the  decoupling  method  of  multivariable  system  design. 

In  the  time  domain,  the  decoupling  problem  has  been  studied 
via  state-space  techniques  by  several  authors  [1-5] .  The  con¬ 
ditions  for  the  existence  of  a  decoupling  system  have  been  devel¬ 
oped  in  these  pioneering  works  [1-5] .  The  state-space  techniques 
are  concerned  with  the  internal  structure  of  the  multivariable 
system.  Thus,  the  limitations  of  the  decoupling  approach  can  be 
derived,  and  a  simple  state  feedback  controller  can  be  designed 
to  achieve  an  optimal  result.  However,  for  a  real  system,  many 
of  the  states  are  not  accessible.  As  a  result,  a  high  degree 
observer  is  often  required  for  practical  application  of  the  state- 
space  technique. 

In  the  frequency  domain,  decoupling  via  the  use  of  transfer 
function  matrices  has  been  investigated  by  several  pioneers  [6-10] 
as  long  as  twenty  years  ago.  Since  then,  practicing  engineers 
have  successfully  extended  the  classical  frequency-domain  approach, 
(for  example,  the  Nyquist  method  [11-13]  and  the  root-locus  method 
[14-15])  from  single-variable  system  to  multi-variable  systems. 

Most  existing  frequency-domain  design  methods  for  multivariable 
systems  either  neglect  the  effects  of  weekly  interacting  . subsystems 
or  completely  destroy  the  coupling  effects  of  the  original  multi- 
variable  system  such  that  a  simple  classical  single-variable  method 
can  be  applied  to  achieve  reasonable  design  goals.  The  approach 
that  removes  the  interactions  of  the  coupled  system  and  designs 
a  controller  for  each  decoupled  system  by  using  a  compensator 
obtained  from  the  plant  inverse  matrix  is  called  the  direct- 
decoupling  method.  This  method  is  straight  forward  but  the 
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following  problems  [16]  arise: 


ik 


a.  The  existence  of  the  plant  inverse  matrix 

b.  The  realizability  of  the  obtained  high-degree  controller 

c.  The  stability  of  the  designed  system  when  unstable 
pole-zero  cancellation  has  been  used 

d.  The  design  procedures  for  a  high-degree  coupled  system 
is  complex. 

Since  many  practical  systems  [12,13,15,17]  are  invertible, 
controllable,  observable,  and  open-loop  stable,  the  direct- 
decoupling  method  is  worthy  of  improvement  and  modification  for 
practical  application.  Recently,  Peczkowski  and  Sain  [17]  have 
improved  the  method  via  a  time-domain  model  reduction  method  and 
have  successfully  applied  it  to  design  the  control  system  for  a 
F-100  engine.  The  direct-decoupling  method  is  modified  and  im¬ 
proved  in  this  paper.  Our  method  uses  the  adjoint  matrix  instead 
of  the  inverse  matrix  of  the  plant  matrix,  and  utilizes  a  fre¬ 
quency-domain  model-reduction  method  for  a  class  of  multivariable 
control  system  design. 

II.  MODEL  REDUCTION  OF  A  MULTIVARIABLE  SYSTEM 

The  usual  procedures  for  designing  high-order  linearized 
multivariable  systems  are  cumbersome  and  computationally  dif¬ 
ficult.  The  controller  obtained  will  usually  be  a  high-degree 
dynamic  controller.  To  overcome  these  difficulties,  a  reduced- 
order  model  of  the  original  system  is  necessary.  Recently, 
various  model  reduction  methods  in  the  frequency  domain  [18-21] 
have  been  proposed  for  determination  of  simplified  models  or  for 
estimation  of  location  of  the  approximate  dominant  poles  and 
zeros  of  the  original  systems.  Since  each  linearized  model  or 
reduced-order  model  is  only  an  approximate  model  of  the  actual 
system,  we  may  use  a  reduced-order  model  as  a  reference  for  the 
control  system  design.  Thus,  the  design  procedure  is  greatly 
simplified,  and  a  satisfactory  lower-degree  controller  can  be 
obtained.  In  this  paper,  various  mixed  methods  will  be  given 
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for  obtaining  the  reduced-order  model  of  the  original  multi- 
variable  control  system.  The  mixed  methods  are  described  in  the 
following  paragraphs. 

Consider  a  typical  transfer  function  matrix  of  a  multivariable 
control  system 


G0(S) 


dc(s) 


(s) 


do(s) 


(1) 


where  each  polynomial  4>  .  .  ( s )  is  the  element  at  the  ith-row  and 

1 • J 

jth-column  in  $  (s)  ,  and  dQ(s)  is  the  least  common  denominator  of 

the  elements  of  $  (s) .  The  transfer  function  g.  .  (s) ,  which  is  an 

1  •  3 

element  in  GQ(s),  can  be  expanded  into  a  continued  fraction  of  the 
Cauer  second  form  [18]  by  using  repeated  long  division  of  the  two 
polynomials  to  determine  the  various  reduced-order  models,  that 
is 


4>,  ^(s)  (b  +b,  s+. .  .+b  sm)  .  . 

g.  .(s)  .  ifj  «  -°  J- - gL-_ig 

1  d  s) 


ao+ais+. . .+ans 


(2a) 


hl  + 


h2  + 


(2b) 


hl  +  hf 


hlh2+s 


1 


(2c) 


hl  + 


h2  + 


h3  +  hf 


h2hsh4+ (h2+h4) s 


^1^2^3^4+  ^1^2+^1^4+^3^4^  s+sJ 


(2d) 
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It  has  been  shown  that  the  first  several  hfs  are  the  dominant 
quotients  in  the  expansion  if  the  steady-state  response  is  of 
interest.  Bosley  et  al.  [22]  have  linked  the  quotients  h^  and 
the  time  moments  (or  the  moments)  of  the  original  system.  The 
reduced-order  models  in  Eq.^2)give  good  approximate  steady-state 
responses.  The  disadvantages  of  the  above  method  is  that  the 
reduced-degree  models  may  not  be  stable  even  when  the  original 
system  is  stable.  To  obtain  the  stable  reduced-degree  model 
of  each  element  in  GQ(s)  and  to  have  the  same  least  common  de¬ 
nominator  polynomial  in  the  reduced-degree  multivariable  system, 
we  may  use  the  following  mixed  methods. 

Let  the  reduced  model  of  the  original  system  in  Eq.  (1)  be 


Gq(s)  = 


do(s> 


$  (S)  = 


*4  (s)  I  * 


where 


*  *i  i(s) 

g.  4  (s)  =  -2^- 


k  k  k 

b  +b.  s+. . . +b^  .s 
o  1  p-l 


p-1 


gi,j(S) 


LH 


(3a) 


Vs> 


k  k  k  p 

a  +a, s+. . . +a  s^ 
o  1  p 


«  1  (3b) 


The  relationships  of  the  dominant  quotients  h.  in  Eq.  (2)  and 

*  *  1  ,  .  . . 
the  coefficients  and  b^  in  Eq.  (3b)  can  be  expressed  in  the 

following  matrix  equation  [23]  by 


[b]  =  [H]  [a] 

where 

•p  *  *  * 

[a]  =  [ , a^ , • . . , ^p_ i^ 

T  *  *  * 

[b]  =  [ bQ ,  b^  , . . . ,  bp_jJ 

(II]  =  [H2]"1[H1] 


(4) 


where  T  designates  transpose,  and 


(H2] 


10  0.00 
0  10.00 

0  0  1.00 

0  0  0  .  0  0 


0  0  0  .  0  hx 


W*! 


10  0.00 
0  10.00 

0  0  1.00 

0  0  0  .  0  0 


0  0  0 


0  h. 


The  h^  are  obtained  from  Eq.  (2)  and  a^  can  be  determined  from  the 

dominant  poles  of  dQ(s)  [21] ,  or  from  the  Routh  table  suggested  by 

Hutton  and  Friedland  [20]  ,  or  the  Routh  table  of  Gustafson  [24]  . 

*  * 

Once  h.  and  a.  are  known,  the  b.  can  be  determined  from  Eq.  (4) . 

1  *  1 

Thus  the  polynomial  j  (s)  can  be  determined,  and  the  equivalent 

dominant  zeros  with  the  preassigned  dominant  poles  obtained  from 

* 

dQ(s)  can  also  be  determined.  Since  each  subsystem  in  the  reduced- 

order  model  has  the  dominant  quotients  h^  (or  equivalent  moments) 

and  the  dominant  poles  of  the  subsystem  in  the  original  system, 

the  reduced-order  model  closely  approximates  the  transient  and 

steady-state  responses  of  the  original  system.  It  also  has  the 

same  tracking  properties  of  the  original  system.  Thus,  the  reduced- 
* 

order  model  G  (s)  is  a  good  approximation  of  the  original  system 
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* 

Gq(s).  Using  G^fs)  as  a  design  reference  model,  the  design  pro¬ 
cess  can  be  greatly  simplified  and  a  satisfactory  low-degree 
compensator  obtained.  This  designed  system  retains  the  equivalent 
transmission  zeros  [25-27]  of  the  original  system.  Our  method 
is  described  in  the  following  section. 


III.  A  MODIFIED  DIRECT-DECOUPLING  METHOD 


Let  the  open-loop  transfer  function  matrix  GQ(s)  of  a  unit- 
matrix  feedback  multivariable  system  be 


Gq(s) 


d0(s) 


$  (S) 


(5) 


where  GQ(s)  e  R(s)mxm  is  a  proper  transfer  function  matrix  and 
dQ(s)  e  R(s)  is  the  least  common  denominator  polynomial  of 
Gq(s)  with  degree  n.  The  expression  $(s)  e  R[s]mxm  is  the  numer¬ 
ator  matrix  polynomial.  Applying  the  first  cascade  precompensator 
(s)  to  the  Gq(s)  yields 


Q1(s)  =  G0(s)K1(s)  =  MsJK^s)  =  diag  {  (6a) 

where 

K^s)  =  adj  $(s)  (6b) 


and 

P (s)  =  det  $ (s) 


(6c) 


Note  that  p(s)  is  a  numerator  polynomial  in  Eq.  (6a).  If  some 
zeros  of  p(s)  are  in  the  right  half-plane  and  all  zeros  of  dQ(s) 
in  the  left  half-plane,  a  practical  compensator  can  be  designed 
without  using  unstable  pole-zero  cancellation.  On  the  other 
hand,  if  the  direct-decoupling  method  is  applied,  the  p(s)  will 
be  a  denominator  polynomial  in  Eq.  (6a)  .  Thus  the  decoupled 
system  is  unstable  and  the  impractical  pole-zero  cancellation 
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approach  would  be  required  to  stabilize  the  system.  Thus,  an 
obvious  advantage  of  our  method  over  the  old  method.  Employing 
the  second  cascade  compensator  K2(s)  to  Q^(s),  we  have  the  diag¬ 
onalized  open-loop  transfer  function  matrix 

f  p(s)k.n.  (s)  | 

Gd(s)  =  QilsjKjts)  .  diag  {  dJTsTdTTs)  )  <’*> 

where 

( kini(s)  | 

K2  —  diag  |  ^ |  ,  i  =  1,2, ...,m  (7b) 

Each  is  an  undetermined  gain  at  the  ith  diagonal  element  of 
Gd(s)  for  the  use  of  the  root-locus  method.  Each  n^ (s)  and  d^(s) 
is  a  scalar  polynomial  to  be  assigned  in  the  design  process.  The 
assignment  of  n^(s)  and  d^(s)  shall  improve  the  performance  of  the 
designed  system  with  the  constraint  that  the  cascade  compensator 
K(s)  -  K1(s)K2(s)  be  realizable.  The  choice  of  n^(s)  and  d^(s) 
is  a  design  freedom  and  experience  is  helpful.  The  total  compen¬ 
sator  becomes 


K(s)  =  K1(s)K2(s) 


kini (s) 

adj  $  (s) diag 


(8) 


Notice  that  G^(s)  in  Eq.  (7)  retains  some  of  the  transmission 
zeros  of  the  original  system.  This  can  be  shown  as  follows: 


Let  Gq(s)  =  $(s)  =  N^sJD^s)"1 

where  NY(s)  e  R(s]mxin  and  DY(s)  e  Rts]™411*  are 
prime  matrix  polynomials.  The  characteristic 
the  zeros  of  det  DY(s)  =  0.  The  transmission 
the  zeros  of  det  NY(s)  =  0. 


(9) 

a  pair  of  relatively 
poles  of  Gq(s)  are 
zeros  of  GQ(s)  are 


Eq.  (9)  can  be  written  as 


Gc(s) 


-  Ms) 

-  do(s) 


Ny (s) adj  Dy(s) 
Mi) 


(10a) 
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where 


A  (s)  =  det  (s) 

*  the  characteristic  polynomial  of  GQ(s) 


=  the  least  common-denominator  polynomial  of  all  minors 
of  GQ(s) .  (10b) 


Taking  the  determinants  of  both  sides  of  Eq.  (10a)  yields 


det 


f*  (SL.  1  = 

Ldo<s>J 


p(s)  det  (s)  det  [ad  j  PY(s)] 


dom(s) 


Am(s) 


AB>~,1(s)det  Ny(s) 
Am(s) 


det  Ny (s) 
A  (s) 


(ID 


Rearranging  Eq.  (11)  gives 


dQ(s) 


A(s) 


det  (s) 


(12) 


Substituting  Eq.  (12)  into  Eq.  (7a)  we  have 

m-1 


Gd(s) 


(k.n.(s)c 
=  drag  J - ^ 


(s) detNfs) 


(s)  A  (s) 
and  the  closed-loop  system  is 

Y  (s)  =  diagj 


(13a) 


kini (s) d0m-1 (s) det  N  (s) 


m-1 


(s)ni(s)det  (s) 


R(s)  (13b) 


(A(s)d.(s)  +  k.dQ 

where  Y(s)  and  R(s)  are  the  output  and  input  vectors,  respectively. 


Since  the  transmission  zeroes  of  the  original  system  GQ(s)  in 
Eq.  (5)  are  the  zeros  of  det  (s)  which  appear  in  both  Gd(s)  in 
Eq.  (13a)  and  the  designed  closed-loop  system  in  Eq.  (13b) ,  the 
designed  system  retains  some  of  the  invariant  transmission  zeros 
when  A(s)  and  det  N^(s)  have  no  common  factors.  In  general,  there 
may  exist  common  factors  between  A(s)  and  det  N^(s).  Therefore, 
Gd(s)  may  retain  only  a  subset  of  the  transmission  zeros.  For  the 


f 


2 

special  case,  when  m=2,  A  (s)  =  dQ  (s) ,  n;.  (s)  =  1,  and  d^s)  = 
dQ(s)  the  compensated  system  Gd(s)  =  (diag{k.det  Ny(s)/det  Dy  (s) }  ) 
in  Eq.  (13a)  has  the  exact  transmission  zeros  as  well  as  the  char¬ 
acteristic  poles  of  the  original  system  G  (s)  . 

When  the  original  system  G  (s)  has  a  high-degree  transfer 

®  it 

function  matrix,  the  reduced-order  model  Gq(s)  in  Eq.  (3a)  can 

be  used.  Thus,  the  designed  controller  K(s)  is  low-degree  and 

the  designed  system  retains  a  subset  of  equivalent  transmission 

zeros  of  the  original  system.  Note  that,  even  if  A (s)  and 

det  N  (s)  in  Eq.  (13a)  have  common  factors,  the  A*(s)  and  det 
*  * 

Ny(s)  obtained  from  the  reduced-degree  model  G*(s)  may  have  no 
common  factor. 


IV.  ILLUSTRATIVE  EXAMPLES 


•Example  1 


To  illustrate  the  design  procedure,  we  use  Mueller's  two- 
shaft  aircraft  gas  turbine  [28]  example.  The  open-loop  transfer 
function  matrix  GQ(s)  of  the  example  unit-matrix  feedback  system 
is 


Go<s>  -  tr^T 


*  (s) 


dQ(s) 


4>u(s) 

<j>2l(s) 


$12(s) 

$22  ^ s\ 


(14) 


where 


dQ(s)  =  s4  +  113.225s3  +  1357.277s2  +  3502.95s  +  2526.85 

4>i;l  (s)  =  14.9s2  +  1506.472s  +  2543.2 

4>12  (s)  =  95150s2  +  1132094.7s  +  1805947 

<|>21  (s)  =  85.2s2  +  8642.888s  +  12268.8 

<J>22  (s)  =  124000s2  +  1492588s  +  2525880 


Notice  that  dQ(s) ,  the  least  common  denominator  polynomial,  is 
the  characteristic  polynomial  of  GQ(s),  or  dQ(s)  =  A(s).  Also, 
this  system  has  no  finite  transmission  zeros.  McMorran  [12] 


designed  a  dynamic  compensator  for  this  system  using  Rosenbrock's 
inverse  Nyquist  array  method  [11] .  His  design  goals  are  a  weakly 
coupled  system  and  a  fast  response.  The  compensator  obtained  was 


Ke(s) 


--26.847s(s+20.6) 

1 

s (s+158 . 5) 

0.018468s  (s+146.3) 


46.272 (s+1. 706) (s+11.6) 
-0.00726  (s+1. 706)  (s+101.4). 


(15) 


Unit-step  responses  are  shown  in  Figure  1. 

The  goals  of  our  design  are: 

(a)  Two  identical  diagonal  subsystems  decoupled  in  the 
closed-loop  system 

(b)  Unity  final  values  of  unit-step  responses 

(c)  Less  than  10  percent  maximum  overshoot 

(d)  The  time  required  for  the  unit-step  response  to  reach 
the  first  peak  of  the  overshoot  tp  is  0.01  sec. 

For  this  low-order  system,  use  of  the  reduced-degree  model  is  not 
necessary.  To  simplify  the  controller  in  Eq.  (8)  let  n^[s)  =  1. 
Then  K(s)  becomes 


K  (s)  =  adj  $>  (s)  K2  (s)  = 


^22  ^  “$12  ^ 
-4>21(s)  4>1X  (s) 


r  k. 


d^s) 


d2(s)J 


(16a) 


The  designed  open-loop  system  is 


Cd(s) 


|-  k1p(s) 
d^7sTd]ls) 


0 


0 

k2p(s) 

dQ(s)d2(s) 


(16b) 


where  p(s)  =  det  $(s)  =  -6259180do(s) ;  kj  and  dA (s)  are  to  be  deter¬ 
mined.  To  circumvent  a  negative  p(s)  in  the  open-loop  transfer 
function,  Gd(s),  we  apply  a  unit-matrix  controller  KQ  to  4>(s) 
before  using  the  controller  K(s)  in  Eq.  (16a) .  Our  objective  is 


To  satisfy  the  first  specification,  we  let  =  k2  =  k  and  d^s)  = 
d2(s)  =  d(s),  thereby  reducing  the  multivariable  design  to  a  scalar 
design  problem.  In  other  words,  we  design  gd (s)  =  6259180k/d (s) , 
a  diagonal  entry  in  Gdn)(s),  to  satisfy  other  assigned  specifications. 
To  meet  the  requirement  of  unit-final  value,  we  choose  gd  (s)  to  be 
"Type  1";  and  to  meet  the  other  two  conditions,  overshoot  and  the  peak 
time,  we  choose  a  second-order  compensator 


gd(s) 


=  6259180k 
s(s+c) 


(17a) 


The  characteristic  polynomial  of  the  designed  closed-loop  system  is 


Ld(s)  =  s2  +  cs  +  6259180k  =  s2  +  2£u>ns  +  un2 


(17b) 


2 

where  c  =  25u)n  and  wn  =  6259180k.  The  parameters  £  (the  damping 
ratio)  and  a>n  (the  undamped  natural  frequency)  are  to  be  determined. 
From  a  design  rule  of  thumb  [29]  we  can  estimate  uin  knowing  tp 

un  =  f"  =  =  300  rad/sec  (17c) 

P 

Also,  another  rule  of  thumb  [29]  gives 


or 


K 


MP 

IT 


Un  0.1  - 

3.14 


(17d) 


Substituting  £  and  un  into  Eq.  (17b)  we  can  solve  for  c  and  k; 
c  =  2£u)n  =  450 

and 


k  «  0.01438  (17e) 

The  characteristic  polynomial  Ad(s)  and  its  poles  are 

Ad(s)  =  s2  +  450s  +  90000 

and 

s1  2  =  -5%  +  junA-Z2  =  -225  +  j!98. 43  (17f) 


The  compensator  in  Eq.  (16e)  becomes 


Vs)  "  i7i+J50T 


‘ -1782. 98s2- 21461. 74s-36319. 33  1368.15s2+16278.3s+25967.5  1 


Ll. 22 508s 2+ 124. 275 s+176. 4116 


The  decoupled  closed-loop  system  is 

90000 


Y(s)  =  diag 


s“  +  450s  +  90000 


R(s) 


-0.214245s-21. 66138 s-36 .5684. 

(18a) 


(18b) 
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The  unit-step  responses  of  the  designed  system  are  shown  in  Fig¬ 
ure  1.  The  final  values  of  the  unit-step  responses  are  unity, 
and  the  maximum  percentage  overshoot  is  1  percent.  Also,  the  peak 
time  tp  is  0.014  sec.  Note  that  the  designed  system  and  the  orig¬ 
inal  system  have  no  finite  transmission  zeros.  Comparing  the  pro¬ 
posed  compensator  Km(s)  in  Eq.  (18a)  and  that  of  McMorran  in  Eq. 
(15)  we  observe  that  both  controllers  are  second  order.  However, 
our  controller  satisfies  more  sophisticated  control  specifications 
than  McMorran* s.  The  unit-step  response  curves  show  that  our 
design  gives  less  overshoot  and  less  oscillations  and  is  completely 
decoupled. 


Example  2. 


To  illustrate  our  design  procedure  using  a  reduced-order  model, 
we  use  a  paper-making  machine  [13]  example.  The  open-loop  transfer 
function  matrix  GQ(s)  of  the  unit-matrix  feedback  system  is 


pn<s)  4>12  (s) 

Go‘s)  *os7*<s>  -  iarW  4  ,  .  .  ,  , 

o  1>  *n(s)  *22(s> 


where 

d-^s)  =  s6  +  34.9798s5  +  565.584s4  +  5016.37s3  +  24517.51s2 

+  55613.33s  +  12868.37 

*11(s)  =  -9.72727s (s2  +  15.39326s  +  112.3596) 

4>i2  (s)  =  173 . 386s  (s2  +  11.44444s  +  55.55556) 

<t»21(s)  =  0 . 204545  (s2  +  15.39326s  +  112.3596) 

4>22  (s)  =  19 . Os  (s2  +  11.44444s  +  55.55556) 


Sinha  and  Rutherford  [13]  have  used  Rosenbrock's  inverse  Nyquist 
array  technique  [11]  to  design  a  controller  for  this  system.  The 
precontroller  to  GQ(s)  is 


Kp(s)  =  J 


m  1  ■ 
s 


-15 (s+1)  52(s+l. 11111) 


3.75 (s+1 . 11111) 
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The  unit-step  responses  of  the  closed-loop  system  using  Kp(s)  are 
shown  in  Figure  2. 

In  this  problem,  our  design  goals  are: 

(a)  The  designed  closed-loop  system  must  be  weakly  decoupled 
and  have  two  approximately  identical  diagonal  subsystems 

(b)  The  final  values  of  the  unit-step  responses  must  be  unity 

(c)  The  tracking  error  should  decay  at  least  as  fast  as  e-*^. 
Gq(s)  is  a  high  degree  transfer-function  matrix.  To  determine  a  low- 
degree  controller  and  to  simplify  the  design  procedures,  a  reduced- 
order  model  will  be  used.  The  continued  fraction  mixed  method  [18] 
and  the  Routh  approximation  method  [20]  may  be  used  to  determine  the 
reduced-order  model  GQ(s).  The  procedures  to  apply  the  mixed  method 
are  as  follows:  First,  the  reciprocal  polynomial  with  respect  to 
d^s)  in  Eq.  (19)  (that  is,  d^fs)  ■  Jd^^))  is  arranged  into  the 
Routh  array  [20]  to  determine  a-^.  The  obtained  is  equal  to 
0.23139.  Thus,  the  equivalent  least  common  denominator  polynomial 
d*(s)  of  the  G*  (s)  is 

d*  (s)  =  sts+a^  =  s(s+0. 23139)  (21) 

Then,  various  first  dominant  quotients  (defined  as  { h ^ ^  ^)  of  each 
entry  in  GQ(s)  can  be  determined  by  performing  the  continued  fraction 
expansion  on  each  ^  (s)  and  dQ(s)  as  shown  in  Eq.  (2).  These  values 
are 

{ hi> 1,1  =  -11.68221 

(hl} 1 , 2  =  i-31079 

{ hi }  2 , 2  =  U-962  (22a) 

and  the  { h^^ } 2 , 1  for  s4>2i(s)  and  dQ(s)  is 

{hl}2  i  =  555.5556  (22b) 

A  *  f  , 

Substituting  the  ai  value  (*=a  )  in  Eq.  (21)  and  each  { h,  }  .  •  in  Eqs. 

X  i  ,  J 

* 

(22)  into  Eq.  (4)  yields  each  reduced-order  polynomial  4» ^  f  j  (s)  .  Thus, 
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Ms) 


(23) 


G0{s)  = 


do<s> 


where 


II 

H* 

{ s) 

4>i2Ms)- 

dO<  s> 

_4>2i*(s) 

*22*<s>- 

,*(s) 

=  s(s+0. 23139) 

o 

11  (S) 

=  -0.01981s 

*21(s) 

=  0.0004165 

12(s) 

=  0.17653s 

f22(s) 

=  0.019344s 

* 

The  unit-step  response  of  GQ(s)  and  GQ(s)  are  shown  in  Figure  3. 

The  approximation  is  excellent.  Using  GQ(s)  as  a  plant  and  closely 
following  the  procedures  used  in  Example  1,  we  have  the  modified 
controller  in  the  form  of  Eqs.  (16)  .  The  precontroller  is 


k1n1  (s)<t>22  (s)  k2n2  (s)  <t>12 (s)" 


Km(s)  =  Koad^  Vs)K2(s)  = 


(s) 

k^i  (s)d>21(s) 

djli) 


d2(s) 

^2n2  ^*11  ^ 


(24) 


d2(s) 


The  designed  open-loop  system  is 

-kip*(s)ni(s) 


Gdm(s)  =idia9  * 
am  1  d0(s)  d.  (s) 


=  diag 


0.0003832 (s+0. 191868) kini (s) 


(s+0. 23139)  di(s) 


for  i  =  1,2  (25) 


where  p*(s)  =  det  $*(s)  =  -0 . 0003832s (s+0 . 191868) 


To  simplify  the  design  procedure,  we  let  n^(s)  =  n2(s)  =  1; 
and  to  meet  the  first  design  goal  we  let  k^  =  k2  =  k  and  d^(s)  = 
d2(s)  =  d(s).  Thus,  we  have  a  single  open-loop  transfer  function 


_  0.0003832 (s+0. 191868) k 
9d(SJ  "  (s+0. 23139)  d(s) 


(26) 
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From  basic  root-locus  theory  (30]  ,  we  observe  that  there  exists  a 
nearby  open-loop  pole,  (pQ  =  -0.23139)  and  zero,  (zQ  =  -0.191868) , 
in  gd(s).  When  gain  k  is  increased,  a  pole  of  the  closed-loop 
system  will  migrate  from  pQ  to  zQ  which  is  a  closed-loop  zero 
also.  Thus,  the  performance  of  the  closed-loop  system  is  heavily 
dependent  on  the  zeros  of  d(s).  To  satisfy  the  third  specification 
we  choose  d(s)  =  s+10  such  that  the  tracking  error  of  the  designed 
system  will  decay  at  least  as  fast  as  e~10t.  The  choice  of  gain  k 
is  a  design  freedom  which  can  ensure  less  influence  of  the  nearby 
pQ  and  zQ.  In  this  example,  we  choose  k  =  40.  The  overall  com¬ 
pensator  in  Eq.  (24)  becomes 


(s) 


1 

s+10 


-0.77376s 
.  0.016667 


7.0612s 

0.7924s 


(27a) 


Finally,  to  achieve  the  second  design  goal  we  add  a  forward-gain 
matrix  H  as  shown  in  Figure  4.  The  H  is 


9d  to) 

0 

-l 

787.786  0 

i+gd(o) 

0 

gd(°) 

1+gd (o)J 

0  787.786 

The  unit-step  responses  of  the  closed-loop  designed  systems  using 
the  compensators  Km(s)  and  H  in  Eqs.  (27)  and  the  reduced-order 
model  G*(s)  in  Eq.  (23)  are  shown  in  Figure  2.  Also,  the  unit- 
step  responses  of  the  closed-loop  systems  using  the  same  compen¬ 
sators  in  Eqs.  (27)  and  the  original  system  GQ(s)  in  Eq.  (19)  are 
shown  in  the  same  figure.  The  response  curves  show  that  the  . 
designed  system  has  less  overshoot,  is  less  oscillatory,  and  is 

nearly  completely  decoupled.  The  discrepancy  between  the  response 

* 

curves  of  the  closed-loop  designed  systems  using  GQ(s)  and  Gq(s) 
appears  in  the  transient  regions.  This  is  the  result  of  the  model 
reduction  method  used.  Both  curves  are  seen  to  be  practically 
decoupled  and  unit-final  values  in  the  steady-state. 


1 
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V.  CONCLUSION 


A  simple  design  technique  has  been  given  for  the  design 
of  multivariable  control  systems.  The  designed  system  retains 
some  or  all  of  the  invariant  transmission  zeros  of  the  original 
system.  When  the  given  multivariable  control  system  has  a  high- 
degree  denominator  and  high-degree  numerator,  various  mixed 
methods  have  been  illustrated  for  determining  reduced-order  models 
By  using  the  obtained  reduced-degree  model  as  a  reference  model, 
the  design  procedures  are  greatly  simplified,  and  a  satisfactory 
low-degree  compensator  has  been  designed.  The  design  technique 
is  simple  and  quite  satisfactory  as  shown  in  the  two  illustrative 
examples.  However,  when  a  multivariable  control  system  has  a 
large  number  of  inputs  and  outputs,  the  degree  of  the  controller 
obtained  by  our  method  may  be  large  despite  using  a  reduced-order 
model.  When  the  resulting  controller  dimension  is  large,  model 
reduction  methods  may  be  utilized  to  reduce  the  controller  to 
an  acceptable  order.  Our  method  is  suitable  for  multivariable 
systems  which  have  high-degree  least  common  denominator  polynomial 
and  low-degree  numerator  polynomials.  Another  advantage  of  our 
modified  method  over  the  original  direct-decoupling  method  is 
that,  for  an  open-loop  stable  multivariable  system,  our  method 
gives  a  simple,  practical  and  realizable  controller  without  using 
the  unstable  pole-zero  cancellation  approach. 
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Figure  1.  Unit  Step  Responses  of  Designed  Systems  in  Example  1 
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ABSTRACT 

An  algebraic  method  is  developed  to  construct  a  new  matrix  Sturm  series 
and  to  establish  a  new  block  canonical  form  of  a  matrix  transfer  function. 
The  matrix  Sturm  sequence  is  then  used  to  determine  the  number  of  real  poles 
of  a  matrix  transfer  function  that  may  not  consist  of  a  pair  of  relatively 
prime  matrix  polynomials,  and  it  can  he  used  to  determine  whether  an  imped¬ 
ance  matrix  ran  be  realized  by  using  KC  elements.  The  block  canonical  form 
is  used  to  construct  a  new  block  state  equation  in  the  block  tridiagonal 
form  and  to  obtain  a  pair  of  relatively  left  prime  matrix  polynomials  of  a 
matrix  transfer  function. 


1.  1 11L  rmliii- t  ion 

The  matrix  transfer  function  T(s)  of  a  multivariable  system  wild  an  m*l 
input  vector  R(s)  and  an  mxl  output  vector  y(s)  is  written  as 

y(s)  =  T(s)  R(s)  (la) 


and 


T(s)  =  1)  (s)  1D2(s) 


(lb) 


where 


Dj  (s) 
D2(s) 


A  s  +A  s  +  .  .  .  +A„s+A  = 

n+1  n  21 

B  s“-1+B  ,s"'2+. . .+B„s+B, 

n  n-1  2  1 


D11sn+D10sn“1+...+D1  s+D, 

11  12  1 ,n  1 ,n+l 


n-1  n-2 ,  ,  _  ,  ^ 

=  D  s  +D  ~s  +...+D  .  s+D 

21  22  2, n-1  2,n 


A.(=U,  .)  and  B.(=I)„  .)  are  mxm  matrix  coefficients, 

i  l.J  i  2 ,  .i 

In  a  single-input  single-output  system,  the  T(s)  in  bq.  (1)  can  be  ob¬ 
served  as  a  scalar  transfer  function  or  a  driving-point  impedance  function. 
The  T(s)  can  be  formulated  into  a  sequence  of  polynomials  (Sturm  sequence) 
to  determine  the  common  polynomial  between  D^(s)  and  D^s)  and  the  number  of 
p-  I  poles  of  T(s)  by  using  Sturm's  theorem.2  When  no  common  polynomial  be- 
U'.  .'ii  l>j(s)  and  1)^  ( « )  exists,  the  T(s)  is  a  function  that  consists  of  two 
rein! i ve 1 y  prime  polynomials,  and  the  dynamic  state  equation  constructed 
using,  Uj(s)  and  1)^ ( s )  is  completely  controllable  and  observable.  This  im¬ 
plies  that  the  dynamic  state  equation  Ls  a  minimal  realization  of  T(s).  An 


excellent  survey  on  the  .applications  of  the  Sturm  theorem  can  be  found  in 


2 


Barnet L  anil  Siljak's  work.^ 

3 

Recently  Bitmoad  and  Anderson  have  developed  the  matrix  Cauchy  theorem 
and  a  matrix  Sturm  series  for  multivariable  systems,  and  they  have  discussed 
their  applications  to  circuit  theory.  In  this  paper,  an  alternate  matrix 
Sturm  series  is  developed  and  a  block  canonical  form  of  a  matrix  transfer 
function  is  also  derived.  Some  properties  and  applications  of  the  developed 
Sturm  series  and  canonical  form  to  the  analyses  of  multivariable  systems  are 
discussed . 

1 1 .  A  Matrix  Sturm  Series 

The  T(s)  in  Kq.  (1)  is  a  real  rational  matrix  with  left  matrix  fraction 
decomposition  D^(s)  If  I)2(s)  is  non  singular,  the  inversion  of  T(s) 

is  expressed  as 

T ( s ) ~ 1  =  D2(s)_1Di(s)  =  sK1+D2(s)"1D3(s)  (2a) 

or 

U  x ( s )  =  D2(s)sK1+D3(s)  (2b) 

wlie  re 

K.  ~  1  1).,  (s)  'i).(k)  as  s  >“> 

Is/  I 

an  in v in  nousingulnr  mil  ri  x  ipiot  ient 

and 


=  Tiie  matrix  polynomial  in  s  witli  degree  (n-1) 


i>2(s)  and  I) 
inversion  of 


D.3(8) 


or 


1)3(8)  = 


where 


“2  *  ft3 
=  an 


and 


I',  (k)  = 
A 


Subst  i  tut  ing 


>:  »3/"J 

j=i  J,j 


The  left  remainder  of  D^(s) 


(s)  have  the  same  degree  of  (n-1).  If  D^(s)  is  nonsingular,  the 
D2(s)  ^D.j(s)  in  Eq.  (2a)  becomes 


L>2  («)  =  »2+D3(s) 


-1 


D4(s) 


(2c) 


D2(s)H2‘-D4(s))l'1 


(2d) 


(s)  * D 2  ( h )  as  s  ’-® 

m*m  nonsingular  matrix  quotient 


Tiie  matrix  polynomial  in  s  with  degree  (n-2) 


Kq.  (2d)  into  Eq.  (2b)  gives 


whe  re 


L)  (s)  =  D./sHsK.+lCS-D.  (s)ll'1 
l  2  L  2  4  2 

=  D2(s)QJ(s)-D4(s)li21 


4 

(2e) 


Q2(s)  =  sK^+ll”1 


In  the  same  fashion,  we  can  use  D^Cs)  with  degree  of  (n-1),  and  D^fs)  witli 
degree  of  (n-2),  to  generate  another  matrix  polynomial  D^(s)  as 


»  (s)  =  04(s)Q3(s)-D6(s)H4 


-1 


(2  f ) 


where 


sK.yt-ll 


-1 

4 


In  general  we  have 


s 


I> 


0  is  an  ni-'m  null  matrix.  The  matrix  coefficients  of  the  matrix  polynomial 

L).(s)  are  tie  lined  as  1).  j  =  1,2,...  .  The  procedure  to  evaluate  the 
t  t  >  J 

matrix  coefficients  D.  .  can  be  easily  accomplished  by  using  the  newly  de- 

i » J 

ve loped  block  Kouth  array  with  block  kouth  algorithm  that  is  different  from 

4 

the  matrix  Kouth  array  with  matrix  Kouth  algorithm  developed  for  expanding 


matrix  continued  fractions. 


The  block  Kouth  array  is  listed  as  follows. 


Ilu-  procedure  to  construct  t  lit*  block  Routh  array  is  described  as  follows. 

Kilter  the  matrix  coef t ic ients  of  D^(s)  and  of  D^fs)  respectively  as  row 
l  and  row  1 ,  and  evaluate  in  succession  by  using  tlie  block  Routh  algoritlim 
(to  be  shown)  to  obtain  row  3  from  rows  l  and  2,  and  row  4  from  rows  2  and 
3.  Then  rows  2  and  4  are  used  to  generate  row  5,  also  the  rows  4  and  5  are 
used  to  generate  row  6.  The  above  processes  to  generate  row  4  from  rows  2 
and  3,  and  row  5  from  rows  2  and  4  are  repeated  and  continued  to  the  last 
row  (2n+l).  The  block  Routh  algorithms  is 


(i)  =  ,J2  I  !>  1 1  *  ranl<  l)r/ 1  =  m 


l>3,.j  =  Yj+rVj+iV  >  l'2, — 11 


-i 


,  =  l)-j||>2l.  rank  =  m 


31 


'^.j  ~  1)2.j+rl>3,  )+lH2’  ]  1.2, . . .  ,n-l 


-1 


(Li)  K  «  .1).  rank  D  ,  “  m  X 

i+ J  i+2 ,1  i,l  i+2 , 1  \ 


lfi  H,j  ''i.j+l  ‘)i+2,j+lKi+l 


\ 

/ 


i  =  1,2,3,. , 
i  =  2,4,6, . 


,  2n-2 


'  i+2  =  "  i+  3,  1  '*  i+2  ,  1  ’  r;mk  Y-3,1  =  m 


i+4,  j  *>i+2,j+l  '  i+3,.j+l"i+2 


(4b) 


0 ,  .  ,  .  ~  0 

■  ’ll  I  ,  ,  I  III 


•r(im  the  block  Routh  array  we  can  construct  a  sequence  of  matrix  polynomial;: 


I).  (s)  ,  i  =  0,1,2, ...  , 


and  we  shall  show  that  the  matrix  sequence 


8 

l!)()(.s)  ,U9(s)  ,D/((s)  ....  ,  U2|i_2  ( «)  }  is  a  matrix  Sturm  series  of  T(s)  . 

i'lie  matrix  series  in  liq .  (3)  can  be  modified  and  expressed  in  the  form 

3 

of  a  matrix  Sturm  series  developed  by  Bitmead  and  Anderson  as: 


U"(s)  =  u’l+2(s)Qm(.s)-D.+4(s),  i  =  0,2,4 . 2n-2  (5a) 

n"  ,..(s)  =  0 
2n+2  m 


where 


Vs)  =  Vs)  =  V's) 


»2(s)  =  0(h) 


%<»>  58  04k(s)|  I!  II  I  \  k  -  -1.2,3, 

i=i 

k  -i 

=  'W>(s)l  11  V,1  ’  k  =  1,2,3, 

j  =  l  J 


(5b) 


and 


Ql(s)  =  Qj(k) 

=  ll2Q-3<«) 


k 

k 

_  1 

'W“'  ■ 

1  II 
j=l 

V,IW“)I" 

)=l 

H4j-2!  ’ 

k  =  1,2,3, 

, 

k 

k- 

1  , 

1  II 
i  i 

V?"!*k-I(,‘)l  11 

11,-1  , 

1  A| 

k  =  2,3,.. 

•)i  +  l(s)  - 

sKi  +  l 

"Cr 

i  =  0,2,4,.., 

,,2n-2 

For  example,  when  i  =  4  in  Fc| .  (5a),  we  have 


9 


04(s) 


l>fi(s)ci5<s)-U8(s) 


(6a) 


From  Eqs.  (5b)  and  (5c)  we  have 


U"(s)  =  D.  (s)H.,1 
4  4  2 

V“>  -  Vs)lC 

"“(«>  •  (66) 


and 


Q”(s)  =  ll/,Q1.(s)II 


-1 


(6c) 


Substituting  Eqs.  (6b)  and  (6c)  into  Kq .  (6a)  and  simplifying  it  yields 


Vs)ll2l  =  |1)6(s)H4ll[H4Q5(s)H21]_U8(s)H61|,21 


(6d) 


or 


U  (s)  =  I)  (s)()r  (s)-l)  (s)II  1 
6  .)  8  6 


(6e) 


(•><>)  is  one  of  the  matrix  equations  in  Kq.  (  i)  .  When  i  =  6,  we  have 


«n) 


Substituting  the  corresponding  D^(s)  and  Q  (s)  in  Kq .  (5)  into  Kq.  (6f)  we 
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have 


Vs)l,4l  =  IU8(s)“61h21,U12,,6Q7(s),141,-D10<s)H81,,41  (6«> 

or 

D6(s)  =  D8(s)Q7(s)-D10(s)H~1  (6h) 

Again,  Eq.  (bli)  is  a  matrix  equation  in  Eq .  (3). 

■>'< 

From  Kqs.  (5)  and  (6)  we  observe  that  each  matrix  polynomial  D^(s)  in 

•k  >V  k 

the  sequence  of  real  matrix  polynomials  {l)_(s)  ,  D_(s),  D,.  0(s)}  in 

U  Z  Zn-Z 

I'.q.  (5)  is  different  from  the  D^(s)  in  {D^  (s)  ,0^  (s)  , .  .  .  ,  ^2n-2  ^  in  O) 

by  various  weighting  constant  matrices  shown  in  Eq .  (5b).  Furthermore ,  botli 

* 

0.  (s)  and  Q.  (s)  are  matrix  polynomials,  and  D.t_(s)  and  l).l0(s)  are  non- 
j  i  i+2  i+2 

*/f 

singular  matrix  polynomials  with  degree  [det  D^+^(s)]  =  degree  [dot  D^+^(s)] 

<  degree  [det  l).,„(s)]  =  degree  [det  D.,„(s)].  Using  bitmead  and  Anderson's 
i.+z  l+Z 

theorem,  we  can  conclude  that  the  sequence  of  real  matrix  polynomials 

k  >'<  * 

I  U(j(s)  ,  U.;  (s)  ,  .  .  .  ,  1)  9(s)}  as  well  as  the  sequence  {U^  (s)  ,  (s)  ,  .  .  .  ,  ^2n-2 

is  a  matrix  Sturm  sequence. 

Ill  '>ist  ri hut  ion  of  Heal  Foots 

It  is  well  known  that  the  scalar  Sturm  sequence  is  often  used  to  deter¬ 
mine  the  number  of  real  roots  of  a  scalar  polynomial  on  the  real  axis  in  the 
i  •  i  *inj » I « •  x  plane.  Hit  mead  and  Anderson  have  extended  l  lie  scalar  Sturm  theorem 
to  the  matrix  Sturm  theorem.  Since  the  matrix  sequence  il)^(s)  ,1)  (s)  , .  .  .  , 

I) ,  ,(s)  [  in  Eq.  (3)  is  a  matrix  Sturm  sequence,  the  matrix  Sturm  theorem 

lu- /. 


can  lie  applied  to  determine  the  number  of  real  poles  of  a  matrix  transfer 


turn:  lion  or  a  rea  I  roots  of  a  matrix  polynomial.  A  simple  criterion  using 
t  lie  block  quotients  K,  and  11^  in  Eq.  (4)  is  developed  to  determine  tlie 
realization  of  an  impedance  matrix  of  RC  networks. b  The  matrix  Cauchy  in¬ 
dex"*  of  a  symmetric  real  raLional  matrix  T(s)  with  T  *(s)  exists  and 
del  T(0)  £  (J  can  he  written  as 


i't(s)  =  Al’[b''(s)_1D"(s)3+Ab[D"(s)  V(s)  ]+.  .  .+AblD*  0(s)~V  (s)  ] 

a  a  0  z  a  2  4  a  2n-2  2n 

=  Ab(l)  (s)"’l)  (s)  l+AblU..(s)~1l)/(s)H:i]+Ab[Hol)/(s)"10,(s)H:1  ] 
a  U  l  a  l  4  l  a  2  4  f>  4 

f  A  b  [  U  1)~  1  (s)  U  (s)  H~ 1  H~  1  |+...  (7; 

a  4  o  o  o  l 


where 


Aa|ni(s)  '  bi+2  (s)  1  =  -21,,(Ti(b))"0(Ti(a))l 


=  Half  the  total  changes  of  signature  of  T.(s)  over 
(u,h).  1  (7b) 

9  +  - 

o(T)  is  defined  as  the  signature  of  a  matrix  T  and  o(T)  =  dim  V  -dim  V 

where  dim  V  is  the  largest  possible  dimension  of  any  subspace  on  which  T 

is  positive  definite  and  dim  V  is  the  largest  dimension  of  any  subspacc  on 

which  T  is  negative  definite. 

Tin*  I  T(s)  in  I'.ci.  (7)  can  be  expressed  in  terms  of  K.  and  li .  obtained 
a  ii 

>v 

in  i.q.  (4)  as  follows.  betting  a  =  -  <>■  and  b  =  ••■and  substituting  IK(s)  in 
I'q.  (A)  into  bq .  ( /)  wc  iiave 

I1’  T  ( s )  =  o(K  )-I+o(ll.,K_)  -  '+a(H  .  K  ll~*)~  *+o  (ll.ll  KJb*)"*+.  .  . 
a  1  IS  4  i  L  lb  7  4 


i  =  l  f  3 , 3, 


«(M.  ) 

i 
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Wile  t'C 


M1  =  K. 


M>  =  "/l 

M4P+1  =  1  "  H4j  IK4S4H-1 (  "i  H43-2I  ’ 

J  1  J 

H«-i  ■ 1  j,  S-j'St-i'j!,  V'1’ 


l  =  1,2,3,. 

i  =  2,3,.. 


(8b) 


when  ,i  =  -  o’  and  b=0,  we  have 


lj|  T(s)  =  i[a(K^)"1+a(K1)"1l  +  ^(H^rVaOl^) 

L'  1 1  '  TU  III.  IX -ii, 

4  5  < 

2n—  1  ...  -i  2n-l  _i 

!  I  )’  (dM?"1)  +  )‘  o(M  )) 

2  i=l,3,5,  1  1=1 ,3,5, 


-1  - 


+  ][a(H4K*ll“1)  i+a(H4K3H21)“1]+... 


(8c) 


who  re 


=  K1  =  *’i  ,n')]  ,n+! 


M1  =  ll2K3 

II 

n„(i)T  no..  ) 

2  4, n-1  2 ,n 

0 

f. 

II 

+ 

* 

1! 

j  = 

U4j  '■K4«.+lt  1  U4j-2 

l  J- 1 

... 

P, 

*v 

1 

11 

11 

i  = 

,  VV.1,!,  "4J 

,-l 


K2  i+l  *’2(1+1),  (n-  i  )  **2  i ,  (n+1-1)  ’ 


t  =  1,2,3,... 
i  =  2,3,... 
i  =  1,2,3,... , n-1 


(8d) 


M.  .ire  shown  in  Mc|.  (Hl>).  When  n=0  and  b= ",  wo  have 


lh  T(s) 


+ 


2 

1 

2 

1 

2 


M^)  ^odc*)'1]  +  -J  lo(H2K3)"1-a(ll2K*)"1l 

la(H4K5H21)"J-o(H4K*ll“1)~1]+... 

2n-l  _  2n- 1  A 

i  }  a(M.  )  -  l  a  (M*. )  ) 

i=l,i.3,  1  i=l,3,5,  1 
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(8e) 


fh  and  arc  shown  in  Eqs.  (8b)  and  (8d). 

The  Cauchy  Indices  1^  T(s)  in  Eq.  (8)  indicate  the  number  by  real  poles 
of  T(s)  in  (a,b).  When  2r  m*m  block  quotients  K.^  and  K.  are  obtained  in  the 
block  Kouth  array  (that  is  constructed  using  an  nth  degree  T(s)  witli  T  ''(s) 
exists  and  det  T(0)  4  0)  and  the  Cauchy  index  in  Eq .  (8c)  equals  to  rm  (where 
r-  n) ,  then  T(s)  is  asymptotically  stable  and  the  multivariable  system  is  an 
aperiodic  system.^  Furthermore,  if  ail  IK^Cs)  ^D.(s)  are  symmetric  real 

rational  matrices  of  an  aperiodic  system,  then  T(s)  is  relizable  as  the  im- 

3  * 

pedance  of  an  m-port  RC  network.  This  implies  that,  when  r  matrices  Q.(s) 

in  Eq.  (5c.)  are  symmetric,  and  the  Cauchy  index  in  Eq .  (8c)  is  equal  to  rm, 

T(s)  can  he  synthesized  using  RC  elements. 

When  the  number  of  distinct  real  roots  of  det  (s)  is  of  interest,  a 

matrix  polynomial  D.;(s)(=  d[D^(s)]/ds)  is  constructed.  By  using  T(s)=D^(s)  ^ 

l)^(s)  and  the  above  procedures,  the  number  of  distinct  real  roots  of  det  L)^(s) 

in  (a,b)  can  he  determined. 


IV  A  lUoi^k  Canonical  Form 

When  the  total  number  (2r)  of  the  block  quotients  K.  and  11.  in  Eq .  (4) 
equals  to  2n,  where  n  is  the  degree  of  b|(s)  in  the  real  rational  matrix 
transfer  function  T(s)(=  l)^  (s)  (s)),  the  left  matrix  fraction  decomposi- 

_  j  ^ 

tion  I)  (s)  l)^ ( w )  in  Eq .  (1)  consists  of  left  coprime  matrix  polynomials, 

which  may  contain  non-symmet ric  matrix  coefficients.  The  T(s)  can  he  for- 
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mul  ail'd  into  a  block  canonical  form  as  follows. 
From  K<| .  (  3)  wo  have 


u2(s)_1n1(s)  =  Q}  (s)-l)2(s)  1 


»4(s)  ' 0  2 ( s )  =  Q3(s)-D4(s)  1 


l>6(s)  lD4(s)  =  Q5(s)-D6(s)-:1 


U2n(s)'  ,)2n-2(s)  =  Q2n-l(s) 


l,2,rtl(l,)'l|>2„<s)  '  "2„ 


1).  (s)H 
4 

U6(s)ll 
Dg  (s)  II 


-1 

2 

-l 

4 

-1 

6 


(9) 


Succ.essivolv  substituting  Fq.  (9)  into  T(s)  yields 

T ( s )  =  1),  (s)_1U2(s)  =  [D2(s)~1D1(s)  j"1 
=  [Q1(s)-D2(s)"1D4(s)H21J"1 

=  [q1(s)-ii)4(s)"1d2(S)  ]-1h21  r1 

=  [Q1(s)-lq3(s)-lU6(s)"1D4(s)  r1!!"1!'1!!"1]-1 

-  iql(H)-ivi3(8)-iQ5(8)-i . . •-iq2n-i(s) ]”1h2.V •  •  ]~]\l ]_1 

(LO) 

wl'  '  I'O 


q. , .  (s)  =  sk  +u 
l+l  l+l  i+2 


i  =  0,2,4,, 


Kq.  (10)  is  a  block  canonical  form  of  T(s).  The  corresponding  block  diagram 


is  shown  in  l'ig.  1  anil  tin.*  block  state  equation  obtained  from  the  block  dia¬ 
gram  can  be  written  us 


z  =  (.;<£  +  hr 

y  =  IV.  (11) 


where 


-(K.,  ,11 

2n-  1 

2n}  a'2n-3H2n-2) 

0 

m 

0 

m 

0 

m 

0 

m 

K2n-I 

_(K2n-3ll2n-2) 

0 

m 

0 

m 

0 

m 

0 

m 

0 

m 

K2n-3 

0 

m 

0 

m 

0 

m 

0 

m 

c;  = 

• 

• 

• 

• 

• 

• 

0 

m 

0 

Ill 

-(K?H 

«7H6)'1 

0 

m 

0 

m 

-  -  -  - 

—  — 

—  —  . 

0 

m 

0 

m 

s1 

1 

i  " 

(w-1 

✓“-N 

1  71 

X 

i> 

w 

1 

1  _°» 

0 

m 

0 

m 

0 

m 

i 

i 

s1  |  - 

(Kj'V‘ 

1  (Ki 

0 

m 

0 

m 

0 

m 

i 

i 

0  1 

m  , 

t 

K-‘ 

i 

i  _(Ki 

i 

Z  -  l/l 

T  T . 

,z2’  ■  ■  ‘  ’  Zti 

kt  =  10 

m 

,0  . 0  ,  1  1 

m  m  m 

• 

f  =  |0 

m 

.0  ,...,0  ,k"‘] 

m  m  1 

The  • 

1'  in  Kq . 

(II)  dec  i  glia  1  es  transpose. 

/. .  a  rc* 

i 

in*  1  ve 

clors,  r  is 

an  iti*l 

i  nput 

ver I  or  and  y 

i  s  an  nr<  1  out  put  vor  lor. 

(i  is  a 

b  1  oek 

L  r id iagona 1 

mat r  ix 

with 

block  element 

s  constructed  by  using  K, 

and  H . 

« 

Fiij.  1.  B1J<‘ 
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An  a 1  tei  nate  block  state  equation  tliat  can  be  directly  written  from 

T(s)  in  Eq .  (I)  with  =  1  is 

11  m 


x  =  Ax  +  Br 


y  = 

Cx, 

x(0)  = 

0 

nmxi 

(12) 

who  re 

0 

m 

0 

m 

0 

m 

!1l,n+l 

D21  " 

X1 

1 

m 

0 

m 

0 

m 

"Dl.n 

U22 

X2 

A  = 

0 

m 

1 

in 

0 

m 

'Ul,n-1 

,  B  = 

1}2  3 

,  x  = 

X3 

0 

m 

0 

m 

I 

m 

-Dl,2 

U2,n 

X 

n 

C  = 

10  . 
m 

0  , 
m 

.  0  , 
m 

lm] 

0nmx  1  an  lln'xl  null  vector  and  is  an  mxm  identity  matrix,  x.  are  mxl 
vectors.  The  block  slate  equation  in  Eq.  (12)  is  an  observable  block  Com¬ 
panion  form.  It  is  interesting  to  notice  that  the  block  linear  transforma¬ 
tion  matrix  between  the  block  coordinates  z  in  Eq.  (11)  and  x  in  Eq.  (12) 
cal,  be  directly  written  from  the  block  Kouth  array  in  Eq .  (4)  as 


x  -  Kz 


03) 
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2n ,  1 

U2n-2,2 

•  U8,n-1 

°6,n-2 

°4,n-l 

U2,n 

ra 

^*2n-2 , 1 

*  D8,n-4 

U6,n-3 

U4,n-2 

U2,n-1 

in 

0 

m 

•  U81 

D62 

»4  3 

U24 

m 

0 

m 

r 

.  0 

m  | 

°61 

1 

M  1 

1 

W  1 

in 

0 

m 

.  o  1 

m 

o  r 

“  1 

U41 

U22 

t 

0 

1 

.  0  1 

0  I 

0  1 

U„, 

m 

m 

m 

m 

1 

m  , 
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Substituting  Kq .  (13)  inLo  F.q.  (11)  and  comparing  the  respective  system  ma¬ 
trix,  input  vector  and  output  vector  of  iiqs.  (11)  and  (12),  we  have 


KCK-1 

(14a) 

KK 

(14b) 

fr"  1 

(14c) 

When  the  total,  number  (2r)  of  the  block  quotients  K_  and  H  is  equal  to  2n, 
the  block  state  equations  in  Kqs.  (11)  and  (12)  are  the  minimal  realizations 
oT  T(s). 

When  2r-'2n,  T(s)  in  Kq  -  (1)  that  does  not  consist  of  coprime  matrix 
po I vnomin 1 s  can  be  written  as: 


T(s)  = 


!>,(«) 


,-l 


.  .+1) 


2, 


=  H:(s)i>t(s)  | 

=  l’1(s)'II>2(s) 


|C(s)l*2(s)  I 


■  l''Llsr+---+|,l,m1'1|l'21sr''+,’22“"2+---+l'2.r1 

(15) 
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determined  l>v  comparing  l he  matrix  coefficients  «f  tlie  following  matrix  eq- 
uat ion : 


D  (s)  =  C(s)P  (s) 


(16c) 


or 


u,  |s"“l+L)  sll‘2+.  .  .+1)  .  |C  .  .  S"-,"I+...+C1]  X 

21  22  2,n  n+l-r  n-r  1 


ll'2iSr-'«'22Sr-2+...+P2ir) 


LL  is  noticed  that  the  block  linear  transformation  matrix  R  in  Fq.  (13)  can 

be  constructed  using  the  block  elements  in  the  new  block  Kouth  array  that  is 

generated  from  T(s)  =  (s)  V..(s)  with  P,,  =  1  but  not  from  those  in  the 

1  2  11m 

block  Kouth  array  generated  from  T(s)  =  D^(s)  '^(s)*  It  is  also  noted  that, 
K  is  an  upper  block  triangular  matrix,  the  inversion  of  R  can  be  obtained  by 
an  iterative  method.  For  example,  the  inversion  of  an  3m*3m  matrix  R^, 
which  is  obtained  by  partitioning  the  R  in  liq .  (13),  is  required.  The  pro¬ 
duct  of  R  anil  K.j  '  is  written  as 


— 

— 

— 

'Vi 

°42 

V) 

“fil 

X 

y 

I 

in 

0 

m 

0 

m 

0 

m 

U22 

0 

m 

z 

= 

0 

m 

I 

m 

0 

m 

0 

m 

0 

m 

,,2» 

i 

c 

3 

0 

m 

0 

m 

0 

m 

i 

m 

x,y  and  v.  are  unknown  block  elements  to  be  determined.  Fxpanding  Fq .  (13) 
and  solving,  the  resulting  matrix  equations  gives 


v  *  -  a  tax! 


21 

(17b) 


l).  z+D  D  1  =0  «•  z  =  -1)“'D  I)  | 

41  22  21  m  41  22  21 


■  V*  ■  -U;i[U42,,;lD22B^i+C23U2i) 


(17c) 

(17d) 


Thus,  we  can  determine  K.^  . 

When  l’^(s)  and  l'^(s)  are  relatively  prime,  the  characteristic  poles  of 
this  multivariable  system  T(s)  are  the  zeros  of  det  P^(s)  but  not  dot  D^(s), 
and  the  transmission  zeros  of  T(s)  are  the  zeros  of  det  P^(s)  but  not 
det  Dy (s) . 

V .  An  Illustrative  Kxampie 

Given:  A  matrix  transfer  function  of  a  multivariable  system  as 


v(s)  =  T (s) R(s) 


(18) 


whe  re 


r ( s )  =  i)  (s)  lo2(s)  =  fC(s)P1(s))'1IC(s)P2(s)] 


with  ii  =  3  and  in  =  2. 

Det  e  rmi in- : 

(i)  The  number  of  real  poles  of  T(s) . 
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(ii)  The  realizability  of  T(s)  usinj>  an  m-port  1<C  network. 

(iii)  The  common  matrix  polynomial  C(s). 

(iv)  The  minimal  realization  of  T(s). 

(v)  A  pair  of  left  coprime  matrix  polynomials  P^(s)  and  P^Cs). 

(vi)  The  characteristic  poles  and  transmission  zeros  of  this  multi- 
variable'  system. 

To  solve  above  problems  we  construct  the  block  Routh  array  as 


5  2 
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Because  D(  |  ~  0^,  tlio  block  Routh  array  terminates  prematurely,  and  we  have 
2r  (r=2<n=4)=4  block  quotients.  From  the  array  we  can  also  determine 


*  - 1 

K  =  0  D  = 

1  0.4 

0.2  \ 

1 

a 

n 

1  23  14  ' 

V  0.2 

0.6  J 

3  42  23 

2  1 
1  3 


(19b) 


Substituting  R.,  R.  and  11.  into  Kq.  (8c)  yields 

ill 


Jlodl"  1)+a(M"  1)+o(m“1)+o(M31)]  =  A 


(19c) 


whe  re 


M,  =  R  j  = 


M1  -  Kl  = 


0.4 

0.2 

0.2 

0.6 

5 

2) 

2 

i  / 

M3  "  H2K3 


M3  "  H2K3 


2  1 

\ 

1  3 

) 

2.5 

-6 

-6 

14.5 

I  i  ,  .  Kq .  (19c)  we  conclude  that  T(s)  has  four  negative  real  poles  in  (-%0). 

To  determine  the  realizability  of  T(s)  using  passive  RC  network,  we 
have  to  test  whether  T(s)  in  Ktj .  (18)  is  a  symmetric  matrix.  This  can  be 

*  ii 

easily  accomplished  by  checking  the  matrices  l)^(s)  and  Q3(s)  in  Kq.  (5c) 
using  R^  and  ii .  in  Kq.  (19a)  as  follows: 


24 


r  sKi+tl2  =  s 


c  h::) 


-i. 


Q.}(s)  =  H2(sK3+H4  )  =  s 


2.5 

-6  14, 


;m:d 


(19d) 


Both  Qj(s)  and  Q^(s)  are  symmetric,  therefore  T(s)  is  symmetric.  Because 
the  number  of  the  Cauchy  index  in  Eq.  (19c)  equals  to  rm  (=4)  and  T(s)  is 
symmetric,  T(s)  can  be  synthesized  using  passive  RC  elements. 

Since  the  block  Routh  array  in  Eq.  (19a)  terminates  prematurely,  we 
can  write  the  C(s)  in  Eq.  (18)  as 


C(s)  =  C5)s+C52  = 


24  22  \  /  54  52 


10  10 


68  64 


1  0 


0  1 


1  i 


s+ 


24  22 


\  2  2  ^10  10 
(19e) 

To  determine  the  minimal  realization  of  T(s),  we  construct  the  block 


state  equation  in  Eq .  (11)  using  the  K.  and  li,  in  Eq .  (19a)  as 


z  =  0z  +  Er  (19  f) 

y  =  Fz 

who  re 


-<w' 

‘W"' 

n  40 
~  V  1  '30 

58}  (  1 

54/  \-2 

1)' 

1 

•C"  QC 

2M-(  1 

10/  V-2 

V_ 

°2 

f  () 

lM 

_ 

l  0 

o  / 

II 

r— i 

1  r—i 

CM 

o 

II 

(° 

('  '2) 

(  ' 

0  \ 

}  0 

0/ 

\-2  5 

'2 

\  0 

1 ) 

—  - 

_ 
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Eq.  (.ISM)  is  tlu.'  minimal  realization  of  T(s). 

The  left  cop rime  matrix  polynomials  P^(s)  and  P  (s)  can  be  obtained 
from  Eq.  (16).  The  C(s)  is  modified  to  ensure  that  P  =  as  follows 


(19g) 


(19h) 


and 


,(s) 


( 1 9  i ) 


Note  that  tile  matrix  coefficients  in  P^(s)  and  P^fs)  are  not  all  symmetric 
but  T(s)  =  P  (s)  *P  (s)  is  a  symmetric  matrix. 

The  characteristic  poles  of  this  multivariable  system  are  the  zeros  of 
dot  P  (s)  =  0,  or 


s,  =  -0.027  19,  s.,  =  -0.127864,  =  -5.88774  and  s.  =  -193.957. 

L  2  3  4 

lli.  transmission  zeros  of  this  multivariable  system  are  the  zeros  of 
del  Pt>(s)  =  0,  or 


Sj  =  -0.10315  and  s2  =  -193.89685  . 
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V I .  Cone  1  us  i  un 

An  algebraic  method  lias  been  developed  for  constructing  a  matrix  Sturm 
series  and  for  establishing  a  block  canonical  form  of  a  matrix  transfer 
function.  A  simple  and  effective  block  Routh  algorithm  has  been  developed 
to  construct  the  block  Routh  array  and  the  block  quotients.  The  block  quo¬ 
tients  have  been  used  to  determine  the  number  of  real  roots  of  a  matrix 
polynomial,  and  to  determine  the  realization  of  a  driving-point  RC  impedance 
matrix.  The  minimal  realizations  of  a  matrix  transfer  function  have  been 
formulated  to  the  block  state  equations  in  the  block  tridiagonal  form  and 
in  the  observable  block  companion  form.  As  a  result,  a  pair  of  left  coprime 
matrix  polynomials  can  be  obtain'd,  and  the  characteristic  poles  as  well  as 
the  transmission  zeros  of  a  multivariable  system  can  be  determined. 
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CHAPTER  V 


CONCLUSION 

The  dominant-data  matching  method  for  analog  pitch  control  system  design 
has  been  successfully  extended  to  design  digital  controller  for  the  semi¬ 
active  terminal  homing  missile  system.  Various  digital  filters  have  been 
designed  and  successfully  tested  in  the  6  Degree-of-Freedom  Terminal  Homing 
Simulation  Program  at  the  MIRADCOM  Laboratories. 

The  developed  direct-decoupling  method  has  been  successfully  applied 
to  design  an  analog  multivariable  gas  turbine  system  and  a  multivariable 
paper  making  machine.  It  is  believed  that  this  method  can  be  further  ex¬ 
tended  to  digital  redesign  of  the  coupled  row  and  yaw  control  system  of  the 
semi-active  terminal  homing  missile  system. 

The  properties  of  the  newly  developed  Sturm  series  and  block  canonical 
form  have  been  discussed.  It  is  believed  that  the  block  canonical  form 
can  be  further  extended  to  synthesizing  a  multi-port  network  function  with¬ 
out  using  integrators. 

Other  new  findings  of  this  research  are  reported  in  the  appendix. 
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Abstract — A  recursive  algorithm  is  developed  for  solving  the  inverse  Laplace  transform,  linear  and 
nonlinear  state  equations  using  block-pulse  functions.  The  relationships  between  the  solution  of  (he 
continuous-time  state  equation  using  block-pulse  functions  and  that  of  the  equivalent  discrete-time  state 
equation  using  trapezoidal  rule  are  investigated.  A  complete  computer  program  is  presented  for  solving  the 
differential  equations  of  linear  and  nonlinear  state  equations  using  block-pulse  functions. 


I  INTRODUCTION 

An  accurate  description  of  a  practical  system  (for  example,  a  semiactive  terminal  homing 
missile  system(IJ)  often  results  in  a  high  order  transfer  function  with  very  large  coefficients 
and/or  a  high  order  linear  and  nonlinear  time-invariant  and/or  time-varying  state  equation  for 
which  the  commonly  used  numerical  integration  methods  (e.g.  the  Runge-Kutta  method |2]) 
may  fail  to  determine  the  time  response.  Recently,  an  alternate  method  (3|  that  uses  the 
block-pulse  functions  has  been  developed  for  solving  the  linear  time-invariant  state  equations. 
In  this  paper,  the  method  due  to  Shieh  [3]  and  others  is  reviewed  and  extended  to  solve  the 
inverse  Laplace  transform,  linear  and  nonlinear  state  equations.  Also,  the  relationships  between 
the  solution  of  the  continuous-time  state  equations  using  block-pulse  functions  and  that  of  the 
equivalent  discrete-time  state  equations  using  the  trapezoidal  rule  [3. 4]  are  further  investigated. 
A  complete  computer  program  based  on  the  proposed  method  is  presented  to  solve  the  inverse 
Laplace  transform,  linear  and  nonlinear  state  equations  using  block-pulse  functions.  Several 
illustrative  numerical  examples  are  included  to  demonstrate  the  superiority  of  the  new  method. 

2.  MAIN  RESULTS 

Consider  a  linear  time-invariant  state  equation 

x(f)  =  Ax[i)+  Bu{t)  (la) 


and 


*(0)  =  xo  (lb) 

where  A  is  an  n  x  n  system  matrix,  B  is  an  n  x  r  input  matrix,  x(t)  is  a  state  vector  of  n 
components,  «(f)  is  a  vector  of  r  input  functions,  and  x(0)  is  the  initial  state  vector.  The 
piecewise-constant  solution  of  (1)  can  be  obtained  by  using  the  block-pulse  functions  <#>,(/)  for 

j  =  1,2 . m.  Each  block-pulse  function  <t>,{t)  is  defined  by  </>;(f)  =  I  for  (j  -  l)T  st  <jT, 

and  <£|(f )  =  0  for  other  cases.  The  term  T  is  a  time  increment  oi  a  sampling  period,  and  m  is  the 
number  of  the  discrete-time  solutions  of  interest.  The  block-pulse  functions  <fy(t)  for  j  = 
1.2. 3. 4  are  shown  in  Fig.  1. 


J 
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The  piecewise-constant  solution  and  the  discrete-time  solution  of  (I)  in  the  interval 
I  j  |)7'  i  <  jT  are  defined  as  the  column  vector  C,  and  the  x*(jT).  respectively.  The  C;  and 
v*(yT)  can  be  determined  from  the  recursive  algorithms  shown  in  the  following  steps: 

Step  I.  Approximate  the  input  vector  u(l)  that  has  r  input  functions  using  the  trapezoidal 

rule.  The  columns  vectors  L,  in  the  r  x  m  matrix  L  =  [JLi,  L; . Lm]  that  is  the  approximate 

input  functions  are 

/.,  -=  *  |«(7'D+  «((/-  nni  for  7=1.2 . m  (2) 

=  average  value  of  //(/lover  the  interval  (/  -  l)T  <  t  <  jl . 

Step  2  Fvaluate  an  n  x  m  matrix  K  =  [K,.  K: . KJ.  The  column  vectors  K,  are 

K,  -  *  BL,  for  j=l.  2 . m.  0) 

Step  V  Determine  an  n  x  m  matrix  D  =  |D,.  D: . Dm}.  The  column  vectors  D,  are: 


D,  =  jR,K, 

(4a) 

0  =  (/„  +  R;)D,  .,+  jR,(K,-  K,  ,)  for  j  =  2.3... 

. .  m 

(4b) 

(4c) 

R  =  R\A 

l„  =  an  «  x  n  identity  matrix. 


Step  4  Obtain  the  n  x  m  required  piecewise-constant  solution  matrix 
C  =  [Ci-  C: . 0,1.  The  column  vectors  C,  are 


C,  =  ^  D,  *  v(0)  <-Sa> 

0  =  0  ,  +  |<D,  ,  +  D,)  for  7  =  2.2 . m.  (5b) 

The  required  piecewise-constant  solution  of  (1)  is 


ittl 
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where  4>  -  (4>i,  <t>-< . 4>m  I'-  The  prime  designates  the  transpose,  and  the  4>,  is  a  block-pulse 

function. 

Step  5.  Determine  the  required  approximate  discrete-time  solution  x*(f)  using  the  reversed 
process  of  the  trapezoidal  rule.  The  required  solution  x*U)  at  discrete-time  t  =  (j  +  1)7"  is 

x*(0  +  l)T)  =  - jc*0T)  +  2C«/ +  1)T)  for  /  =  0.1.2 . m-l 


x*(0)  =  x(0)  and  C((j  +  \)T)  =  C,. ,.  (6) 

The  expression  x*{jT)  is  the  approximate  discrete-time  solution  of  the  x(f)  in  (I).  The 
accuracy  of  the  approximation  depends  heavily  upon  the  chosen  sampling  period  T.  A  complete 
computer  program,  based  on  the  algorithms  in  (2)— (6)  is  presented  in  this  paper  to  obtain  the 
solution  x*(jT)  in  (6). 

Because  the  trapezoidal  rule  is  applied  to  approximate  the  input  function  in  (2).  it  is 
interesting  to  investigate  the  relationships  between  the  solution  of  the  continuous-time  state 
equation  using  block-pulse  functions  and  that  of  an  equivalent  discrete-time  state  equation 
using  the  trapezoidal  rule  [4.  5]. 

When  «(f)  =  0  in  (I),  Shieh  et  a/.(3]  have  shown  that  the  equivalent  discrete-time  state 
equation  of  the  continuous-time  state  equation  in  (1)  is 

x*((j  +  l)T)  =  Gx*(jT)  (7a) 

x*(0)  =  xn 


G  =  /„  +  /?3=/„  +  (/„-|at)  'at  (7b) 

=  (/n  +  ^7')(/„-^r)  =(/„-^A7-)  '(/n  +  ^Ar).  (7c) 

The  discrete-time  solution  of  (7)  is 

x*{jT)  =  G'x(0)  for  /  =  0.1.2...  (8) 

In  this  paper,  the  derivation  is  extended  to  a  more  general  case  (Shieh  et  o/.[3]),  that  is.  u(t)  *  0. 
To  simplify  the  expression,  the  T  in  the  following  discrete-time  state  equations  and  difference 
equations  is  dropped.  Also,  the  vectors  L,  are  expressed  by  L(j).  K,  by  K(j),  D,  by  D(j)  and  C, 
by  C(j). 

When  «(f)  *  0.  (2)  to  (5)  can  be  expressed  by  a  set  of  difference  equations 

D(l)  =  y  R,K(\)  =  j  R,Ax(0)  +  y  R,BL(\)  (9a) 


D(j  +  !)  =  (/„  +  R2)D(j)  +  Y  +  I)  -  Uj)]  for  /=  1.2 . m 


C(l)  =  i(2/„  +  /?j)x(0)  +  ^/?,BL(l) 


Oj  +  I)  =  Cij)  +  y[D(/+l)  +£>(/)]  for  /=1.2....i 
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Substituting  (9)  into  (10)  and  using  (6)  will  yield  the  required  discrete-time  solution  x*(j),  or 
x*(J)  =  4>*U)x(0)+  <*>*(/ -  i)RiBL(i) 

=  4>*(i)x(0)  +  'Z<t>*(i-i-\)Hu*(i)  for  /=  I.2....  (I la) 

i= 0 

where  <t>*(j)  =  The  transition  matrix  of  the  discrete-time  system 

-[('■-H'H'17')]' 

R,=  r(/„-^7-)  ’ 

H  =  RtB 

«*(i)  =  i[«(/+l)+u(t)l  (lib) 

T  =  sampling  period. 

From  (I la)  the  approximate  discrete-time  state  equation  can  be  written  for  the  continuous-time 
state  equation  in  (I)  as 

**(/+!)  =  Gx*(j)+  Hu*(j)  (12a) 

x*(0) =  x(0)  (12b) 

where 

G  =  (/„-^r)  '(/„  +  ^r) 

H  =  \at )  B 

u*(j)  =  |  («(/+  !)+«(/)]•  (12c) 

It  is  believed  that  the  modeling  of  the  discrete-time  state  equation  in  (12)  from  the 
continuous-time  state  equation  using  block-pulse  function  is  new.  If  the  Z  transformation  is 
performed  on  both  (1)  and  (12).  we  have  the  respective  functions  as 

zX(z)-zX(0)=  AX(z)+BU(z)  (13) 

and 


zX*(z)  -  zX*( 0)  =  GX*(z)  +  HU*{z) 


(14) 


where 


l/*(z)  =  J[zl/(z)+l/(z)]. 

Substituting  the  G  and  H  in  (12)  into  (14)  and  simplifying  yields 

f  (TTT)  X*{Z)  ~  f  (771)  (/n  ~  \  AT)X*M  =  ***<*)  +  BUiz)  05) 


■ 


where  X*(0)  =  X(0). 
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Comparing  (13)  and  (15)  and  assuming  that  X*(z)  is  the  approximate  function  of  XU),  we 
have 

Z(i(01  =  f(FTT)X*(z)"f(7TT)(/"-^i47')X(0)  <16) 

Equation  (16)  is  the  approximate  numerical  differentiator  that  is  often  used  to  determine  the 
inverse  Laplace  transform  of  a  continuous-time  state  equation [4, 5).  Thus,  the  solution  of  a 
linear  time  invariant  state  equation  can  be  obtained  from  the  recursive  algorithms  in  (2)-(6).  or 
from  the  matrix  equations  in  (11)  For  linear  and  nonlinear  time-varying  systems,  the  frozen¬ 
time  and  frozen-state  approach  [6]  is  applied  to  solve  the  linear  and  nonlinear  problems  using 
block-pulse  functions.  In  other  words,  when  an  independent  variable  t  and  the  time  dependent 
state  variables  x^t)  appear  in  the  system  matrix  in  (I)  at  stage  j.  the  time  variable  t  is 
considered  as  a  frozen  time  by  letting  t  =  jT.  The  state  variables  considered  as  frozen  states 
Xj(jT)  in  the  time  intervals  jT  1  )T.  Substituting  these  constants  t  =  jT  and  xtU)  = 

xt(jT)  into  the  system  matrix  in  (I)  and  using  x(jT)  as  the  initial  vector  yields  the  time-invariant 
state  equation  in  the  time  interval  jT  ■<  t  s(j  +  \)T.  Thus,  the  proposed  method  can  be  applied 
to  evaluate  the  solution  x(t)  at  f  =  (/+l)T.  Repeating  the  processes  we  have  the  required 
discrete-time  solutions  for  the  linear  and  nonlinear  state  equations.  A  complete  computer 
program  based  upon  the  above  approach  is  presented  in  this  paper  for  solving  the  inverse 
Laplace  transform,  linear  and  nonlinear  state  equations. 


3  FORMULATIONS  OF  STATE  EQUATIONS  AND  ILLUSTRATIVE 
EXAMPLES 

Consider  that  the  impulse  response  of  the  following  rational  function  is  required: 

E(r). 

U{s)  r"  +  Oij"'*  +  a:* •  +  on 

Since  the  input  function  (/($)=  1,  the  required  impulse  response  is  the  inverse  Laplace 
transform  of  Y(s).  Also,  the  impulse  function  is  a  delta  function,  it  cannot  be  realized  because 
of  its  infinite  amplitude  at  t  =  0.  Therefore,  it  is  convenient  to  convert  (17)  into  a  zero-input 
state  equation  with  initial  conditions  as 


jf  (/ )  =  Axit)  +  Bui  l) 

(18a) 

y(t)=  Cx(t) 

(18b) 

x(0)  =  xo 

(18c) 

where 


‘  0  1 

0 

0  - 

'x,(r)i 

~0~ 

0  0 

1 

0 

.  *(/)  = 

•MO 

B  = 

0 

a«  ~  Qn-t 

-  a,_ 

.0 

C=/„. 

The  output  matrix  C  is  chosen  as  an  identity  matrix  so  that  the  output  functions  are  equal  to  the 
state  variables.  The  initial  vector  jr(0)  can  be  evaluated  from  the  following  matrix  equation 
(2,7]: 


x(0)  =  D'b 


(19a) 
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I 


or 


'  v(0) 

*,(0f 

~1 

0 

0 

•  0 

Jt‘(0) 

Jt?<0) 

«i 

1 

0 

0 

v'(0) 

- 

X\(0) 

— 

fl; 

a  i 

1 

0 

x1”  "(0) 

x„(0) 

a„  i 

a„  2 

Qn  f 

l_ 

(I9b) 


b  is  a  vector  constructed  from  the  coefficients  of  the  numerator  polynomial  in  (17).  A  recursive 
algorithm|2|  has  been  proposed  to  determine  the  initial  conditions  without  finding  the  inversion 
of  the  square  matrix  D.  An  alternate  method  is  proposed  in  this  paper  to  determine  the  D' 1 
and  the  required  initial  vector. 

To  determine  the  D  1  indirectly,  we  construct  the  following  matrix  equation: 


z  =  K(-  a) 


(20a) 


or 


*1 

I 

0 

0 

0 

0 

o' 

”-fll~ 

r? 

1 

0 

0 

•  0 

0 

-  O: 

:i 

= 

1 

0 

■  0 

0 

-  fll 

, 

-4 

r  t 

sj 

7  ( 

1 

0 

0 

-  04 

I 

*n  2 

’ft  ^ 

2n  -4 

•  2| 

I  _ 

(20b) 


The  vector  (  -  a)  is  constructed  from  the  coefficients  of  the  denominator  polynomial  in  (17)  with 
negative  signs.  The  matrix  K  is  a  lower  triangular  matrix  with  each  diagonal  entry  assigned  as 
unity,  and  other  entries  2,  are  determined  from  the  vector  z.  In  other  words,  from  the 
multiplication  of  the  first  row  vector  in  K  and  the  vector  (-  a)  we  have  the  numerical  value  2,. 
Then,  we  immediately  substitute  the  z,  into  the  lower  diagonal  entries  and  solve  for  z?.  and  so 
on  The  general  algorithm  is 

:o=  I 

z,  =  -  £  z,  ta,  ,.i  for  j  =  1.2 . n.  (20c) 

I  I 

The  matrix  K  is  the  inversion  of  the  matrix  D  in  (19).  Thus  the  required  initial  vector  x(0)  can 
be  determined  in  ( 19).  If  the  input  vector  «(/)  in  (17)  can  be  expressed  by  analytical  functions  or 
a  set  of  finite  values  at  sampling  points.  (17)  can  be  converted  to  a  zero  initial-state  time- 
invariant  state  equation  as 


x(t)=  Ax(t)  + Bu(t)  (21a) 

y(f)  =  Cx(t)  (21b) 

x(0)  =  0  (21c) 

where  A  is  shown  in  (18);  B  =  [0, 0 . 0,1]’;  C  =  [b„.  b„_i . bi,  M‘<  and  x(0)  = 

(0.0 . 0.01'. 

A  practical  system  is  the  transfer  function  of  the  pitch  control  system  of  a  semi-active  terminal 
homing  missile  system(l]  which  is  shown  as  an  illustrative  example  as  follows: 
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where 


a,=  1.9235540  x  10' 
«;.  =  9.3162391  x  10' 
a,  =  2.9769507  x  10* 
u4  =  6.2316753  x  10'° 
a,  =  9.3603299  x  10'-’ 
<i,  =  9.7499233  x  1014 
a7  =  6.6673970  x  lO16 
a,  =  2.4204054  x  10'* 
u9  =  2.91 19206  x  10'* 
<i,o  =  2.4190474  x  JO1'4 
u,,  =  8.8021585  x  I01* 


=  0 
6,  =  0 
=  0 
6>  =  0 
64  =  0 

6s  -  1 .4945233  >  10" 
6„  =  2.5633964  x  |014 
67  =  5.0172120  x  lO16 
6,  =  2.926344.3  x  10'* 
K  =  4.6100047  x  lO'4 
6,„  =  8.8021585  x  H)1*. 


It  is  desired  to  find  the  step  response  Uis)  =  (1  Is). 

Equation  (22a)  can  be  formulated  either  in  the  form  on  (18)  or  that  of  (21).  Attempts  to  solve 
this  problem  by  the  Runge-Kutta  method (2)  were  unsuccessful  even  though  the  time  increment 
was  chosen  as  small  as  10  4  sec.  This  is  because  the  practical  system  consists  of  large 
coefficients  in  the  transfer  function.  This  normally  results  from  large  poles,  for  example,  in  (22a> 
there  exists  a  small  a,  (which  is  the  sum  of  all  poles)  and  a  large  a ,,  (which  is  the  product  of  all 
poles).  This  difficulty  is  overcome  by  the  proposed  method.  Using  the  proposed  computer 
program  with  time  increment  DT  =  0.2  sec  yields  the  u;ut  step  response  curve  shown  in  f  ig.  2 
For  comparing  the  results  of  the  proposed  method  and  the  Runge-Kutta  method  we  apply  both 
methods  to  the  reduced  third-order  model  of  the  original  llth-order  system  in  (21a)  (using  the 
method  of  Shieh  and  Chen|8.9])  to  evaluate  the  unit-step  responses: 

Y*U)  0.6920s2  +  19,4692s  +  3.7376 

Uis)  ?To.9488r  +  lOlbbTv  +  3.7376 '  ' 

The  response  curves  are  shown  in  Fig.  2.  From  these  results,  we  observe  that  the  proposed 
method  is  superior  to  the  Runge-Kutta  method  if  the  system  consists  of  extremely  large  or 
small  coefficients  or  the  response  curve  of  the  system  has  many  stiff  slopes. 

When  a  linear  or  nonlinear  time-varying  equation  is  given  and  the  numerical  solution  is 
required,  the  given  equation  can  be  converted  into  a  state  equation  in  (1)  with  lime-varying  and 


Fig.  2  Time  responses  of  the  systems  described  in  eqns  (22;;)  and  (22b). 
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nonlinear  entries  in  the  system  matrix.  The  proposed  method,  along  with  the  frozen-time  and 
frozen-state  method,  can  be  applied  to  determine  the  solution.  To  illustrate  the  procedure  we 
use  the  following  examples: 

Given  a  nonlinear  equation 

^  -  /(/)[1  -  z(t)2]  ^  +  Kz(t)  =  <?«(/)  (23a) 

z(0)=  a  i  and  z(0)  =  <*2 


where  /(f)  is  a  time-varying  function,  and  K,  Q  and  a,  are  constants.  Equation  (23a)  can  be 
converted  into  a  state  equation  by  defining  the  state  variables  *,(/)  and  x2(t)  as 

x,(f)  =  z(f) 

x2(t)  =  z(t).  (23b) 

The  corresponding  state  equation  is 


pMOlf  0 

U(/)J  L-K 


/(OH  -  x.2(Ol]  [jrz(o]  +  [q]  "(/) 

KH:;] 


(23c) 

(23d) 


When  f(t)  =  1,  K  -  I  and  0=0,  the  nonlinear  equation  is  the  Van  der  Pol  equation[10].  If  the 
output  functions  y,(t)  and  y 2(<)  are  assigned  as  x,(t)  and  x2(0.  then  the  state  equation  is 


BKH-t  ,-^KMSH 

fy.(01  ri  01fx,(r)l 

U«)J  10  iJlx2(/)J 

rx,(0)]  fa,i 

UWj'UJ 


The  proposed  method  that  uses  the  block-pulse  functions  and  the  frozen-time  and  frozen-state 
approach  can  be  used  to  solve  (24)  for  determining  the  trajectories  x,(t)  and  x2(t).  In  other 
words,  substituting  x(0)  in  (24c)  into  the  system  matrix  in  (24a)  results  in  a  time-invariant  state 
equation  in  the  form  of  (I).  Thus,  the  developed  recursive  algorithm  in  (2)-(6)  can  be  applied  to 
evaluate  x*(T)  that  is  the  required  discrete-time  solution  x(t)  at  f  =  T.  Then,  using  the  x*(T) 
obtained  as  the  new  initial  vector  in  (24c)  and  again  substituting  the  x*(T)  into  the  system  matrix 
in  (24a)  to  obtain  the  new  time-invariant  system  matrix  for  evaluating  the  new  solution  x*(T) 
that  is  the  required  solution  x(t)  at  t  =  27,  and  so  on.  The  trajectories  of  the  nonlinear  equation 
as  a  result  of  different  sets  of  initial  conditions  are  shown  in  Fig.  3. 


4.  CONCLUSION 

The  recursive  algorithm  for  solving  linear  time-invariant  state  equations  using  block-pulse 
functions  has  been  extended  for  solving  the  inverse  Laplace  transform,  linear  and  nonlinear 
time-invariant  and  time-varying  state  equations.  The  relationships  between  the  solution  of  the 
continuous-time  state  equations  using  block-pulse  functions  and  that  of  the  equivalent  discrete- 
time  state  equations  using  trapezoidal  rule  have  been  investigated.  It  is  shown  that  the 
discrete-time  solutions  of  both  methods  are  identical.  An  approximate  numerical  differentiator 
has  also  been  derived.  A  complete  computer  program,  based  on  the  derived  recursive  algorithm 
using  block-pulse  functions  and  the  frozen-time  and  frozen-state  approach,  has  been  written  for 
solving  the  inverse  Laplace  transform,  linear  and  nonlinear  state  equations.  It  has  been  shown 
that  the  proposed  method  is  superior  to  the  Runge-Kutta  method  if  the  system  consists  of  stiff 
functions. 
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APPENDIX 

This  program  is  used  to  solve  the  inverse  Laplace  transform,  linear  and  nonlinear,  time-invariant  and  time-varying  state 
equations.  The  details  to  prepare  the  input  cards  can  be  illustrated  by  the  following  examples. 

Example  I.  For  the  following  transfer  function: 


.  ,  s!  +  8j  +  I2 
y(s>  s^nPTiosTo 


(Al) 


The  discrete-time  responses  y(t)  at  (  =  #T  for  T  =  0.25  sec  and  /  =  0, 1 . 4  are  required. 

The  input  nomenclature  follows: 

The  first  data  card: 

KS— Type  of  problems  to  be  solved.  When  KS  =  1,  it  is  the  inverse  Laplace  transform  problem.  For  this  example,  KS  =  I. 
N— Degree  of  transfer  function.  For  this  example,  it  is  3. 

MT— Number  of  discretr  lime  solutions  required.  In  this  case,  it  is  5. 

DT— Time  increment,  o.  npling  period.  For  this  example,  it  is  0.25  sec. 

The  second  data  card: 

AA0,  AA(I) . AA(N>— Coefficients  of  the  denominator.  For  this  example,  they  are  1. 3,  -  10.0. 

BB(!) . BB(N)— Coefficients  of  the  numerator.  In  this  case,  they  are  1.8, 12. 

If  the  degree  of  the  numerator  is  less  than  (N  -  I),  zero  coefficients  are  assigned  in  the  numerator.  The  output  data  of 
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this  program  arc: 

KS...I  N...3  MT...5  DT. ..  0.250000 


DENOMINATOR  COEFFICIENT  AA0.AAII). .  .AA(N)  ARE 
O.IQOOOOOOE  01  0.30000000E  01  -O.IOOOOOOOE  02  0. 

NUMERATOR  COEFFICIENT  BB( I). .  BB(N)  ARE 
O.IOOOOOOOE  01  0.80000000E  01  O.I2000000E  02 
INITIAL  CONDITIONS  X(l). .  .XIN)  ARE 
O.IOOOOOOOE  01  0.50000000E  01  0.70000000E  01 

THE  REQUIRED  SOLUTION 


T 

Yl 

Y2 

Y3 

0. 

O.IOOOOE  01 

0.50000E  01 

0.70000E  01 

0.25000E  00 

0.25897  E  01 

0.77179E  01 

0.I4744E  02 

0.50000E  00 

0.5I446E  01 

0.I272IE  02 

0.25283E  02 

0.75000E  00 

0.938IOE  01 

0.21 169E  02 

0.42302E  02 

O.IOOOOE  01 

0.I6436E  02 

0.35275E  02 

0.7054 IE  02 

Note  that  VI  =  y(r)  and  V2  =  >(t). 

Example  2.  Consider  the  following  state  equation: 


x(t)  =  Mt)  +  Bu(l) 

y(l)  =  Cx(l)  (A2) 

x(0)  =  x0 


where 


A  =  the  system  matrix  =  f  *  B  =  the  input  matrix  = 

C  -  the  output  matrix  =  ^  ®  j;  a(t)  =  the  input  vector  = 
x(0)  =  the  initial  vector  =  j^J  J;  y(r)  =  the  output  vector® 


The  input  functions  «■(()  and  aj(t)  are  unit-step  functions.  The  responses  y(() at  f  =  jT  for  T  =  0.25  sec  and;  =0, 1 . 4 are 

required. 

The  input  nomenclature  follows: 

The  first  data  card: 

KS— Type  of  problems  to  be  solved.  When  KS  =  2,  it  is  the  problem  of  solving  linear  time-invariant  state  equations;  when 
KS  =  3,  for  solving  linear  time-varying  and  nonlinear  state  equations.  In  this  example  KS  =  2. 

N— Order  of  the  state  equation.  For  this  case,  it  is  2. 

MT— Number  of  discrete-time  solution  required.  In  this  example,  it  is  5. 

DT— Time  increment  or  sampling  period.  For  this  example,  it  is  0.2S  sec. 

The  second  data  card: 

KU— Type  of  input  functions.  If  KU  =  I.  the  input  functions  a,(f)  are  continuous-time  functions.  All  uf(r)  can  be  inserted 
in  the  main  program  using  /  as  an  independent  variable.  If  KU  =  2.  the  input  functions  are  in  the  form  of 
discrete-time  input  data.  For  illustration,  in  this  example.  KU  =  2. 

NU— Number  of  the  input  functions.  For  this  example,  it  is  2. 

NP— Number  of  the  output  functions.  In  this  case,  it  is  2. 

The  third  data  card: 

A(l.  I ).  A(  1 . 2) . A<  I ,  N )— The  entries  of  the  first  row  vector  in  the  system  matrix  A.  For  this  example,  they  are  1.2. 

The  fourth  data  card: 

A(2, 1 ),  A(2. 2) . A(2,  N )— The  entries  of  the  second  row  vector  in  the  same  matrix.  In  this  case,  they  are  3,  -  4. 

The  fifth  data  card: 

B(l,  I).  B(l,2) . B(l,  NU) — ' The  entries  of  the  first  row  vector  in  the  input  matrix  B.  For  this  example,  they  are  2.0. 

The  sixth  data  card: 

B(2, 1),  B(2,2) . B(2.  NU)— The  entries  of  the  second  row  vector  in  the  same  matrix.  In  this  case,  they  are  1.1. 

The  seventh  data  card: 

C(l,  1). C(l. 2) . C(1 .  N ) — ' The  entries  of  the  first  row  vector  in  the  output  matrix  C.  In  this  example,  they  are  1.0. 

The  eighth  data  card: 

C(2. 1).  C(2, 2) . C(2.  N)— The  entries  of  the  second  row  vector  in  the  same  matrix.  For  this  example,  they  are  0. 1. 

The  ninth  data  card: 

x(l) . x(N>— The  initial  conditions.  For  this  example,  they  are  1. 1. 

The  tenth  data  card: 

all.  I), all, 2) . all,  MT)— The  discrete-time  data  of  the  first  input  function  u,(i)  evaluated  at  r  =  ;T  for  T  =  0.25  sec 

and ;  =  0, 1, _ (MT  -  I).  In  this  example,  «,(()  =  I:  therefore,  the  discrete-time  input  data  are  1, 1. 1, 1. 1. 

The  eleventh  data  card: 

a(2, 1).  a(2.2) . a(2,  MT)— The  discrete-time  data  of  the  second  input  function  u2(()  evaluated  at  t  =  ft  for  T  =  0.25 

second  and  ;  =  0, 1 . (MT  - 1).  For  this  example,  «,(()  =  I;  the  discrete-time  input  data  are  !„l,  1, 1, 1. 
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Since  the  Ui(r)  and  u,( t )  are  unit-step  functions  (continuous-time  functions),  we  can  choose  K U  =  I  lf  we  choose  KU  =  I. 
we  do  not  need  the  10th  and  llth  data  cards.  However,  the  statements 


UT(I)  =  I. 
UT(2)  =  I. 


should  be  inserted  in  the  main  program. 

The  output  data  of  this  program  are: 

KS...2  N.  .  .2  MT  .5  DT. .  .  0.250000 
KU.  .2  NU...2  NP...2 

SYSTEM  MATRIX 
O.IOOOOOOOE  01  0.20000000E  01 

0.30000000E  01  -0.40000000E  01 

INPUT  MATRIX 
0.20000000E  01  0. 

O.IOOOOOOOE  01  O.IOOOOOOOE  01 

OUTPUT  MATRIX 
O.IOOOOOOOE  01  0. 

0.  O.IOOOOOOOE  01 

INITIAL  CONDITIONS  X(l) . . .  X(N)  ARE 
O.IOOOOOOOE  01  O.IOOOOOOOE  01 

DISCRETE  INPUT  DATA 
1.000  1000  1.000  1.000  1.000 

1.000  1.000  1.000  1.000  1.000 


THE  REQUIRED  SOLUTION 


T 

Yl 

Y2  Y3 

0. 

O.IOOOOE  01 

O.IOOOOE  01 

0.25000E  00 

0.25897E  01 

0.I564IE  01 

0.50000E  00 

0.51446E  01 

0.27883E  01 

0.75000E  00 

0.938 10E  01 

0.48942E  01 

0.10000E  01 

0.I6436E  02 

0.84191  E  01 

Example  3.  Given 

a  nonlinear  state 

equation  in  (24)  of  the  main  paper 

GKH-:  .-iJBM" 

r>,(ol  ri  oirx.wi  r*,«»]  ui 

ly2(/)J  10  lJU(i)J'  U(0)J  LojJ 

where  i/(r)  =  0. 

The  procedures  to  prepare  the  input  data  cards  are  the  same  as  those  shown  in  Example  2.  except  the  following: 

(i)  KS  =  3. 

(ii)  If  KU  =  I  is  used,  the  statement  of  the  input  function  UT(I)  =  0  is  inserted  in  the  main  program  as  shown  in  the  list 
of  this  program. 

(iii)  Any  entries  that  consist  of  nonlinear  or  lime-varying  terms,  or  both,  in  the  system  matrix  are  first  assumed  to  be 
zeros  in  the  input  cards,  then  the  corresponding  terms  are  recovered  by  the  exact  terms  by  substituting  the  state  variables 
Xjlr)  by  x(j)  and  r  by  T  in  the  main  program.  In  this  example,  the  system  matrix  in  (A3)  is  first  assumed  as 

i]  (A4> 

In  other  words,  the  entry  A(2, 2)  is  a  nonlinear  term  which  is  assumed  to  be  zero  in  (A4).  The  term  is  recovered  in  a 
statement  A(2, 2)  =  I  -x(l)»x(l)  in  the  main  program  as  shown  in  the  list  of  this  program.  The  outputs  of  this  example  are 
plotted  in  Fig.  3  of  the  main  paper.  The  complete  computer  program  for  solving  the  inverse  Laplace  transform,  linear  and 
nonlinear  state  equations  follows. 
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Computer  program 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

5  00 

SOI 

sin 

600 

601 

A.r> 

605 

604 

605 
60* 
607 
60S 
603 
610 
611 

61? 

6?0 


4  PROGRAM  FOR  SOLVING  INVERSE  LAPLACE  TRANSFORM 

LINEAR  STATE  EQUATIONS  AND  NONLINEAR  STATE  EQUATIONS  USING 

PLOCR-PULSf  FUNCTIONS. 

K$»1  INVERSE  LAPLACE  TRANSFORM 

<S*2  LINEAR  TIME-INVARIANT  STATE  EQUATIONS 

<$*?  NONLINEAR  AND  LINEAR  TIME-VARYING  STATE  EQUATIONS 

WHEN  KS*1/ 

Y ( S ) *  <nq< 1 ) *S**<  N-1 > ♦.. .♦ft*MN) )/ <AAO*S**N*AA (1 >«S**  <N- 1 >♦... 

♦  A  A  (  N  )  ) 

WHEN  <S *2/ 

0  X ( T ) /0T«AX<T)4PU(T) 

Y <T)»CX  <T> 

X(0)«  INITIAL  VECTOR 
WHEN  <5*5/ 

THE  TIME-VARYING  AND  NONLINEAR  ENTRIES  IN  THE  MATRIX  '  A  *  ARE 
ASSUMED  TO  OP  ZERO  IN  THE  INPUT  DATA  CARDS/THEN  THE  CORRESPONDING 
ENTRIES  ARE  RECOVERED  BY  THE  EXACT  TERMS  IN  THE  MAIN  PROGRAM 
USING  STATE  VARIABLES  X <  1  >/.../ X ( N )  AND  TIME  VARIABLE  '1'  . 

N*  DEGREE  OF  THE  TRANSFFR  FUNCTION 
*  ORDER  OF  THF  STATE  EQUATION 
MT*  NO.  OF  THE  KNOVN  DISCRETE-TIME  IN»UT  OATA 
*N0.  OF  THE  OUTPUT  DATA  REQUIREO 
OT*  TIME  INCREMENT  OR  THE  SAMPLING  PERIOD 

WHEN  K U *  1 / 

THE  IN°*JT  FUNCTIONS  UT<J>  ARE  C  ON  T  !  NUOU  S  -  T  I  ME  FUNCTIONS. 

WRITE  ALL  THE  L'T(J)  IN  The  MAIN  PROGRAM  USING  THE  INDEPENDENT 
VARIABLE  #  T  * . 

WHEN  <’J*?/ 

THE  INPUT  FUNCTIONS  are  0  I S C» E T F- T I *E  DATA. 


N  I)*  NO. 
NP»  NO, 


of  the  input  functions 

OF  THE  OUTPUT  FUNCTIONS 


F0RMAT(5lS/f 10.6) 

FORMAT  (<5E1*.R>) 

FORMAT ( (AFIO.X)) 

FOR»»AT</?X/SM«(S,../I2/4X#4HN.../!?/CX/5HMT,,./I5,4X/Sh0T.#./F10.6) 


,  A  A  (  N  i 


ARF  •  ) 


«0RMATW/4X/  ♦DENOMINATOR  COEFFICIENT  AAO/AAd). 

FORMAT  I  f?X/6F 16. 8) > 

FORMAT ( / / 4  X / •NUMERATOR  COEFFICIENT  QM < 1 >•••** (N ) 

FORMATC  *6X/'  INITIAL  CONDITIONS  X(1).,.X(N)  ARE’) 

FORMAT (  I  >X/*F 16, « ) ) 

FORMAT(/?X/5H<U.../I2#AX/5hNU.../I2/AX#SHNP. . .  /  !2) 

FORMAT ( //4X, ' SYS TEM  MATRIX') 
f 0RMAT(//6X/’ tN»UT  MATRIX’) 

F  0RM4T  (  / /6X  /  '  OUT  otjT  MATRIX*) 

FORMAT ( //6X/ ' DISCRETE  INPUT  DATA*) 

FORMAT (/ /6X/ * THP  REQUIRED  SOLUTION  •  / 9 X / •  T * / 1 2 X / • Y 1 • 

STTX/,r5,/1TX/,r4,/11X/*Y5t/11X^,v6,/11X/,v7') 

FORMAT ( TX/QE 1  A. S ) 
f  0  R  M  A  T  (  (  ?  X  /  R  p  1  n  ,  5  )  ) 

dimension  aa(20>/5B<20)/A(20/20)/A1(20*20)/XN<20)#x(20)#xG(?0)f 
srl  (??#  100) /0<  20/ ?0)/C  <  20/ 20  )#'.'(  20#  inn)  /  ML  (20/ 100)  /  AG  (  a'' WUT(?0)/ 

«X<1 (70) #  XI <  ?0) ^RT (20/20) /R2 (20/ 20)/ XD( 20/ 100) /XD1 (?0)/xr( 2  0/ 10P)/ 

txO?(?0) ,xo 3(70)/ XC  <20 #110)/ XT (20/100)^7 (?0/ 100)  /R5( 70/21 ) 


ARF  *  ) 


•/I  1 X/ •  Y?' 


c 

c 

moo  read  (S/500/End*i  )  <S/N/MT/nr 

WRITE (6/600)KS/N/MT/DT 
msMT- 1 

IF(*S,Ne.1)  GO  TO  1 00  T 
WRHE(6/601  ) 

pEA0(5/SO1)  » A0# ( a A(L ) /L*1 /N) 

WRITE (6/*0?)  A  AO, ( AA(L) #L*1 #N) 

WRITE(6/603> 

PEAD(S/S01 ) (B°(L)/L*1 #N> 
WRITE(6/602) (°n(L)#L*1/N) 
no  101  l ■ 1 / n 

P9(L)*BB(L)/A*0 

101  AA(L) *- AA(L) /»A0 
Nl *N-1 

DO  102  L *  1 /N1 
L 1 *L ♦  1 

DO  103  J  * 1 / N 
105  A(L/J)*0. 

102  A(L#Ll)*1. 

DO  134  J*1/N 
L  *N*1 - J 

104  A(N/J)*AA(L) 

DO  110  L  * 1 / N 
DO  120  J  * 1 / N 
if (L/J) *0. 


120 
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m  ftlCL#L)»1, 

N  ?  *  N  -  1 
PO  U1 
<  *0. 

no  160  J  s 1 # L 

^(^0  s*ami,,j).aa(j>»$ 

L 1 *t ♦ 1 

00  1  7rt  JJ*Ll#N 

J  Jt»J J-l 

WO  A1(JJ,JjL)«S 

ICO  CONTI NUf 

no  180  1*1, N 

*N<1  )  »0% 

00  180  J*1#N 

*N<L>*Al(l.,j)*nq(j)*XN<l) 

» (L)*xvrL  > 

180  *G<t,)**U> 

W«! T£ (6,604) 

wPtTE(6,605)  <*<LWL*1,N) 

no  ion  l*i,n 

no  100  j  a  1 #  v 
101  *L(L»J>*Oa 
\P*N 

no  70? 

no  ?ni  j«i#v 

?''o  c  ( l /L  )  * i  • 

r,o  to  ?ooo 

1001  RfAP(5,$00>  8U,NU#NP 

w«trc<o#mO)  « u*nu,np 
w«  l  TF (6,607) 
nO  ? 1 1  l ■  1  , M 

»SAM5,S01>  !»(UJ)/JsM) 

?10  ui»1TE<6#002><A(l*J>#J*1#N> 

■a nr  te(6,oop> 
no  7 1 1  l*1  ,n 

REAP<5,5ni><q(l,j>,j*l,Nll> 

?n  wRITF(6,602H«<1.,.>>,J*1,NU> 

an \  re (4, 619) 

nO  ?  1  ?  L  a  1  ,  N  *» 
read(s,soi  Hr(uj)#j*i/N) 

71?  uRITE<n,50?)<Ca#J>/J*1*N) 

wri te  co,ooo 

R€Ap(S,S01 )|I(J)/Ji1/N) 

n o  ?H  c a  1  »  *1 

IN( l) »X ( l 1 

?  U  X  r,  (  L  )  *  x  (  L  ) 

If<* U.E0.7)  oo  TO  251 
T*n. 

no  ?61  J  s  1  ,  M  T 

no  761  LsI/«'J 

c 

c  •••  IF  X  *J*  1  ,  INSERT  Hi  C  ONT  l  NOQUS- T  1 ME  INPUT  FUNCTIONS 

C  UT<n,...  ,UT<N’J)  H  F  R  €  AND  USE  *T*  AS  AN  INDEPENDENT  V»PM9lE.*»* 

C 

'IT  <  U*0. 

C 

c 

c 

2f>  1  nU.#J)«UT(l> 

^6^  r*7+or 

G o  T->  ?S2 

7S1  a*ITF(6iM0> 

DO  715  l»1/NU 

«EA©<5#510><»iU,J>/J*1/*1T> 

715  ««rffM/1?0)(U(L/J)/J*1#8T) 

25?  no  215  l*i»NU 
DO  215  J«1,* 

J  1  »J  *1 

215  UL(L/J»«(Ua*J)*'J(L/Jl))/2, 

no  7D0  x  L  *  1 ,  N 

DO  700 

*•0. 

no  730  K<*T,Nff  . 

5*S»9<«L#*K)*UL(K'(#KJ) 

700  «L(XL/K J)*S 

!M»S«E9.?)  GO  TO  2000 

x*0 

T«0. 

5000  r  OMT I  N'i f 

C 

C  IF  <S*A,  REPLACE  T  WE  ENTRIES  THAT  ARE  NONLINEAR  AND  f  OP  T  I  *8  F 

C  VARYING  TFR^S  IV  THE  MATRIX  »A*  R  Y  THE  EXACT  TERMS  USING  THE  STATE 

C  *M),  ,XCN>  AND  THE  T  I  *E  VARIANCE  'T',*** 

C 

A ( 2/2)«1,'X<1  >*X  (1  > 

c 

C 

c 

T*T  *DT 

<**  *1 
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CALL  *MLTP(N,N,4,1,XG/*G> 

DO  ?1i>  l*t,N 
DO  ?1  7  J  *  1  #  N 
41  <L#  l> «-A(t # J > / >. 

<k  1<L/l  >  *1.  /PT»«1  (L/l  > 

(Ul  D1NV*  <A1#R5#N) 

*0  ?05  L*1,* 

*0  7Q5  J*1/N 

|f  <*X.nf.5>  GO  TO  4000 
00  ? ! 4  L*  t ,H 

XKl(l)*<AG<L>*0L(L/O>/?. 

T4LL  «MlT*(N#N/RT#T/XK1/X1  > 

00  >10  L  *  1 ,H 
*M<L»**»<L>**6<L> 

*(L>«-XG(L)>?.»X05<L> 

*  r,  (  l  )  *  x  c  i  > 

*C (L#f 1 *X0X(L  > 

I  f  (X.EO.N)  oo  TO  5000 

00  TO  5000 

nO  >10  <L*1/*J 

»>n  710  *J«1*N 

5*0. 

r»0  710  <<*1,N 

^Ki^muioo^tK/cj) 

»?(*L#< J)*S 
DO  210  L»1,N 
00  ? 50  J»1#N 

XK <L/ J ) **G(L>  *BL <L/J ) 

DO  ?  5 1  J«1#N 

If  (J.GT.l)  GO  TO  501 

00  ?52  L*1/N 

X  X  1  (l)**K(L#J  >  /"T 

CALL  1ALTB(N#M/R1 »1*XK1;X01 > 

00  ?5>  L*1#N 
X  0  (  L  /  J > -X01 <L  ) 

GO  TO  251 
J2*J-1 

no  254  L*1,N 

XKl(L>aXK(l,J>/DT-XIC<L/J2WDT 

CALL  *ALT0<n,N#P 1#1 /XKl /X02) 

CALL  MMTP(N,N,P2,1,XD1,X05> 

00  25  s  L«1#N 

x0<L/J)=X0(L#J2)>x05(L)^XO?(L) 

X  0 1 <L>*XO'l,j  > 

com  I  Nil  E 

00  ?A0  JS1,M 

IF  (J.GT.l)  fin  TO  502 
00  241  L*1,N 

XC<l*J>*XML/J>*l'T/2.*XG<L> 

GO  TO  >40 
J  2*J-1 

00  24?  L*t,N 

Xf(L/J)*XC<L#j2)>(X0(L#J2)4X0(L#J))*OT/?. 

CONTINUE 
DO  4,00  La1/M 
X  T (  L  /  1  >  SXN(L) 

DO  401  L  * 1  *  N 
00  401  J>?,¥T 
J?*J-1 

XT(L#J)**XT(L#J>)>2.*XC (  L  /  J  ?> 

00  7?0  *l»1/NP 
no  ■’20  KJ«1/*T 
<=0. 

oo  7?n  Kxai/N 
S*S*C<<L/<K>*XT(*X/KJ> 

1 ( <l*K J ) *5 
tfRTTE  <  6#6  1  1  ) 
r*0  40?  J*1#«T 
T« ( J- 1 ) -OT 

WBITF(6#61?>  T#<r<l/J)#L*1/NP> 

GO  TO  T^OO 
STOP 
F  NO 

SUBROUTINE  NALTP<N/N/ A#X,0/C> 

DIMENSION  A(?0/?0>/B(20>/C(20> 

l*«k 

00  10  L*1/N 

s*o. 

00  10  J  ■  1  /  N 
S*$»ACL/J>*B( J> 

C(L)*S 

RETURN 

€N0 


SUBROUTINE  DINVN  <A,D,N> 

0  I  Nf  NS  I  ON  A <  ?P/20  >  /t> (  20/21 ) 

00  10  !*1/N 
00  10  J  *  1  /  N 
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10  at  t.J>*»< I. J> 

*  1 

J  J»0 

do  ?n 

do  ?o  j*i ** 

XP  *><j/«)«0. 

MKe*)*’. 

J J»KtfO 
LL*  J  J 
K<a«KO 

40  I  f  (AOS <D(JJ/*K> > -1 .E-4) SO/ 50/60 

$0  J  J  * J  J ♦ 1 

RO  TO  40 

60  if  tu.-JJ wo,«n/*o 

>o  *»o  00 

t>TFMO«0(LL/WM) 

OTLL/'*  *•)»»(  JJ/**> 

60  6  (  J  J  z**  )  *0  T€**  p 

*P  (IJV'IKM) 

*0  130  lJ*1z" 

J  ***1 -L J 

1«3  0<</J)*P<«zJ>/OIV 

DO  110  I  * 1 # N 
FAC  *0  (I/O 
DP  110  L J  *  1 /* 

J**+1 -LJ 

I*  <  I-<)1?0/1 10/1 ?0 
1 70  D< f/J)«D(t /J)-FAC*0(</J) 

110  C0NMTJE 

DO  130  J»1/N 
1*0  M  J/*' 

?n  '■ONTIN'JC 

<ttrtt7V 
F.*4D 
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Synthesis  of  optimal  block  controllers  for  multivariable 
control  systems  and  its  inverse  optimal-control  problem 

Y .J.  Wai,  M.Sc..  and  Prof.  L.S.  Shiah,  M.Sc.,  Ph.O. 


Indexing  terms:  Multivariable  control  systems.  Control-system  synthesis,  Optimal  control 


Abstract 

A  new  method  is  presented  to  synthesise  optimal  block  controllers  for  a  class  of  multivariable  control  systems 
represented  by  the  block  companion  form.  The  reverse  process  of  obtaining  the  optimal  block  controller  is  used  to 
determine  the  block -weighting  matrices  of  the  quadratic  performance  index  from  prescribed  control  specifications. 


1  Introduction 

Tlte  accurate  description  of  linear  time-invariant  systems  in 
the  time  domain  may  result  in  m  r/th-degree  coupled  differential 
equations,  or  an  nth-degree  matrix  differential  equation  with  m  x  m 
matrix  coefficients'  as 


n  *  1 

£  A^'X  =  U 

iM 

n 

(la) 

y  =  I  CyD'-'x 

1*1 

(lb) 

D*"1  *(0)  =  «„  /  =  1,2,.. 

(lc) 

where  y  is  an  m  x  1  output  vector,  u  is  an  m  x  1  input  vector  and  x  is 
an  m  x  1  state  vector.  A,  and  C,  are  m  x  m  matrix  coefficients,  and 
the  differential  operator  D  =  dldt.  When  each  initial  vector  a,  is  an 
m  x  1  null  vector,  the  corresponding  frequency-domain  representation 
of  eqn.  1  is  an  nth-degree  matrix  transfer  function  written  as 

Y(s)  =  T(s)V(s)  (2a) 


C  =  fC,  C,  ...  CJ  (Ad) 

The  block  elements  A,,  Om.  lm  and  C,  are  m  x  m  constant  matrices. 
m  x  m  null  matrix,  m  x  m  identity  matrix  and  m  x  m  constant 
matrices,  respectively.  The  vector  X  consists  of  n  blocks  ( Xy, 

i=l,2 . n)  and  each  m  x  1  block  X,  consists  of  m  state  variables. 

In  this  paper,  we  define  the  vector  X  as  a  block  vector.  Because  the 
state  equation  in  eqn.  4  is  formulated  in  the  phase-variable  block 
form,  the  X  is  defined  as  a  vector  in  the  phase-variable  block  co¬ 
ordinate.  As  a  result,  the  Y(0)  is  an  initial  block  vector.  From  a 
conventional  viewpoint,  the  same  vector  X  is  viewed  as  a  vector  with 
nm  state  variables  in  a  general  co-ordinate.  Therefore,  the  same  state 
equation  in  eqn.  4  is  viewed  as  a  state  equation  in  a  general  co¬ 
ordinate.  In  this  paper,  all  the  derivations  are  based  on  the  state 
equation  in  the  phase-variable  block  co-ordinate  rather  than  a  general 
co-ordinate. 

Tlte  objectives  of  this  paper  are  described  as  follows: 

(a)  Obtain  the  optimal  block-control  law  u  =  —  R~'BrPX  =  —  KX 
(where  the  feedback-gain  matrix  K  =  R'  BrP  consists  of  m  x  m 
block  elements  K,,  n)  to  minimise  the  quadratic  per¬ 

formance  index 


where  Y(s)  and  U(s)  are  the  m  x  1  output  vector  and  the  m  x  1  input 
vector,  respectively,  and  the  matrix  transfer  function  T(s)  is 

T(s )  =  Nr(s)D;'(s)  =  Dj'(s)N,(s)  (2b) 

The  matrix  polynomials  Dr(s )  and  Nr(s)  with  appropriate  size  are 
right  coprime,  D,(s)  and  N, (s)  left  coprime.  Let  us  define 

Dr(s)  =  lmsn+Ansn-'  +...  +  /M  +  A,  (3) 

Nr(l)  =  O"-'  +  Cn-yS"-1  +  .  .  .  +  C2S  +  C, 

where  A,  and  C,  are  m  x  m  constant  matrices.  The  corresponding  first- 
degree  state  equation  in  the  controllable  phase-variable  block  form  or 
in  the  controllable  block  companion  form  is 

X  =  AX  +  Bu  (Aa) 


•/  =  -(  \XrQX  +  uTRu]dt  (5n) 

2  -o 

for  the  dynamic  system  formulated  in  the  phase-variable  block 
co-ordinate  in  eqn.  4.  Tlte  T  designates  transpose,  the  weighting 
matrix  R  is  an  assigned  m  x  m  positive-definite  matrix,  and  the 
block-weighting  matrix  Q  is  an  assigned  nm  x  nm  nonnegative 
definite-symmetric  matrix  with  m  x  nt  block  elements  Q,j  =  Qj,. 
or 


0.1 

0.2  • 

•  o.„ 

Qu 

Q22  • 

02„ 

<5b) 


y  =  CX:x(0)  =  X„  (A b) 

where 


om 

om 

om 

om 

X, 

A  = 

om 

om 

lm 

om 

.  B  = 

om 

,  x  --- 

X2 

-Ay 

~~  a2 

-A, 

~*n 

lm 

X„ 

(Ac) 
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Qn  1  On  2  -  -  -  Onn 


The  nm  x  nm  matrix  P  is  the  positive-definite  solution  of  the 
steady -state  Riccati  equation2 

PA  +  A  TP  +  Q  -  PBR  'BtP  =  Onm  (5c) 

Tlte  same  P  can  be  also  solved  from  the  following  canonical  form:2 


X 

A 

—  BR~'Br 

X 

G 

-Q 

—  At 

G 

G<«)  =  «(»)  =  Onm*,.X(0)  =  Af0  (5  d) 


nmn _ innno/tvtfuuQ  +  mint-  stiitt _ 


Expanding  eqn.  1 


It  is  noted  that,  if  the  pair  [A,  B]  is  controllable  and  the  pair 
M.  L\  is  observable  (where  Q  =  LLr),  then  the  closed-loop 
system  is  not  only  optimal  but  stable. 

(b)  Determine  the  block-weighting  matrices  Q  and  R  of  the  quad¬ 
ratic  performance  index  in  eqn.  5 a  if  the  optimal  block  controller 
K  is  assigned  or  if  the  closed-loop  poles  (or  the  equivalent  control 
specifications3)  of  the  optimal  controlled  system  are  prescribed. 


0i  —  On  —  Qu  F  Aj  RA  i 

D2  ~  O12  —  Qji  =  O12  +  AlRA2  —Q21  ~  a]  RA 1 


D2n.,  =  R  (10c) 


2  Linear  optimal-block -regulator  problem 


Taking  the  Laplace  transform  of  eqn.  10a  and  neglecting  the  initial 
conditions  we  have  the  matrix  polynomial  0(s): 


In  the  conventional  synthesis  of  the  linear-regulator  problem, 
the  state  equation  in  eqn.  4  is  viewed  as  a  state  equation  in  a  general 
co-ordinate.  An  optimal  control  law  is  then  derived  by  solving  eqns. 
Sc  or  5 d.  In  this  paper,  the  state  equation  in  eqn.  4  is  considered  as  a 
state  equation  in  the  phase-variable  block  co-ordinate.  The  optimal- 
block-control  law  is  derived  as  follows. 

Expanding  eqn.  4  and  adding  a  trivial  identity  yields 

X,  =  AT, 

*.  =  X2 

X,  =  x3  =  x2 

XM  =  xn  =  -/I,  AT,  —A}Xj  — . . .  —A„Xn  +u  ( ba ) 

Rewriting  the  last  equation  in  eqn.  6 a  gives 

u  =  A,  AT,  +A1Xl+...  +  An*r')  +  X\n)  (6  b) 

Substituting  eqn.  6  into  eqn.  5a,  we  have  an  alternate  form  of  the  cost 
function  as 

F(X,  u)  =  F{X , ,  AT, _ X{,n>)  =  F(X')  =  }  x'tQ'x'  (7) 


where 


Qn 

O12  - 

■  ■  oI» 

A^R 

X, 

Q*t 

022  ■ 

Q*n 

aJr 

x, 

0*  = 

Qnl 

On2  - 

Onn 

at„r 

II 

10 
• 

* 

• 

11 

A}"-) 

RA, 

ra2  . 

..  RAn 

R 

Ar}n> 

Qtj  -  Qij  +  Aj RAj  =  QjJ 

The  (n  +  l)m  x  (n  +  l)m  constant  matrix  Q*  is  a  block  weighting 
matrix  with  mxm  block  elements.  Applying  the  gradient  matrix 
operations4  to  the  quadratic  cost  function  in  eqn.  7  yields 

Fx,  =  Vm  Om  ...  Om\QX' 

J'Fx,  =  lOm  lm  O„]0*Af* 


dn 

jp;FX}")  =  [om  om 

...  /m]e*Ar*<", 

(8) 

Substituting  eqn.  8  into  the  following  Euler’s  equation5 

„  d  d2  „ 

(9) 

Fx - Fx  4-  — 5 ~FX  — . . . 

'  dt  *'  dt 2  ' 

+  (-Dn^M">  =  OmX, 

we  have 

D,Af,  +  D2X,  +DjAf{J> +  .. 

where 

.  +  DJn+1Ar<2">  =  omXt 

(10a) 

[O. 

D2 

0, 

■  02n*|] 

=  [/m 

-lm 

fm  ■■■ 

(-1) 

"In, 

Qn 

Qn 

0*3 

..  OTn 

a]r 

0m  . 

om 

Om 

Om 

om 

Q21 

Qn 

•  ■  Os.n-l 

Ql,n 

aIr  . 

■  Om 

Om 

Om 

om 

om 

Qn 

•  •  Qi,n-a 

Oa.n-I 

Ql.n  ■ 

■  om 

Om 

Om 

om 

om 

om 

■  On,! 

On,  2 

On. 3  ■ 

•  Qn,n 

AlR 

Om 

0m 

0m 

om 

■  0m 

RA, 

ra2  . 

•  RAn 

-,RA„ 

R 

(10*) 


0(s)*.(s)  =  [0JnMi2"  +D2ns2"-'  +  ... 

+  0jj+0,U,(s)  =  omX.  (11) 

where  D2kt,  =  02kt2 ,  k  =  0,  1 .  n  and  D2k=—D2k, 

k  =  1 , 2 . n.  It  is  well  known  that  the  poles  of  the  state  equation 

in  eqn.  5 d  are  symmetrically  distributed  about  the  origin  in  the 
s-plane,  so  are  the  roots  of  the  determinant  of  the  matrix  polynomial 
0(s)  in  eqn.  11.  Performing  the  spectral  factorisation6,7  of  the  matrix 
polynomial  0(s )  results  in  a  stable  matrix  polynomial  A(r)  and  an 
unstable  matrix  polynomial  A  (— s),  i.e. 

D{s)  =  FTL(-sflU)F  (12) 

where 

R  =  FtF  =  D2n„ 
and 

A(s)  =  /ms"  +  E„sn-'  +  . . .  +  JT2s  +  f, 

The  required  optimal-block-control  law  is  then  obtained  from  eqns.  bb 
and  1 2  as 

u  =  [X,  (13) 

where 

A,  =  Ai-Ei,  /  =  1.2 . n 


When  the  given  system  is  not  in  a  phase-variable  block  form,  a  newly 
developed  algorithm  shown  in  Appendix  8  can  be  applied  to  obtain  a 
block  linear  transformation  that  transforms  a  class  of  state  equations 
in  a  general  co-ordinate  into  the  phase-variable  block  co-ordinate. 
Thus  the  proposed  method  can  be  applied  to  determine  the  optimal 
block  controller. 


3  Inverse  optimal  control  problem 

Given  a  set  of  prescribed  closed-loop  poles,  or  equivalent 
control  specifications.3  we  wish  to  determine  the  weighting  matrices 
Q  and  R  of  the  quadratic  performance  index  in  eqn.  5a  by  which 
the  controlled  feedback  system  has  prescribed  closed-loop  poles  and 
the  feedback-control  law  is  optimal.  This  is  an  inverse  optimal-control 
problem.  Kalman*  initiated  the  inverse  problem  for  a  linear  time- 
invariant  single-input  system.  Chang.9  Tyler  and  Tuteur'0  have 
studied  the  problem  via  the  root-locus  method,  while  Molinari," 
and  Anderson  and  Shannon12  have  investigated  the  problem  for  a 
multivariable  system.  All  the  developed  methods  are  based  on  the 
system  equation  formulated  in  a  general  co-ordinate  rather  than  in  a 
phase-variable  block  co-ordinate.  Since  the  multivariable  dynamic 
system  is  formulated  in  a  matrix  differential  equation,  it  is  more 
natural  to  investigate  the  problem  in  the  phase-variable  block  co¬ 
ordinate  than  that  in  the  general  co-ordinate. 

It  is  well  known  that  a  feedback-gain  matrix  can  always  be  ob¬ 
tained  to  give  a  system  with  prescribed  closed-loop  poles  if  a  system 
is  controllable.  However,  the  feedback  controller  may  not  be  optimal. 
In  this  paper  we  determine  the  block-weighting  matrices  Q  and  R  of 
the  quadratic  performance  index  by  which  the  feedback  controller 
not  only  provides  the  controlled  system  with  prescribed  closed-loop 
poles  but  also  performs  optimally.  The  steps  involved  are  described  as 
follows: 

Step  I 

Define  a  characteristic  matrix  polynomial  A(s)  of  the  desired  closed- 
loop  system  whose  matrix  coefficients  consist  of  some  unknown 
parameters  (for  example,  the  damping  ratio  |  and  the  undamped 
natural  angular  frequency  oj„  etc.)  to  be  adjusted.  TheA(s)  is 

A(s)  =  /„sn  +  Ens"-'  +  ...  +  E2s  +  Et  (Ho) 

If  the  desired  characteristic  polynomial  of  the  closed-loop  sv<em  is 
[rf(s) lm  =  U"  +</„*"■'  +  . .  ,+d2s  +rf,)m  (14b) 
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“"whore  </(s)  is  a  polynomial  whose  coefficients  consist  of  adjustable 
parameters.  The  characteristic  matrix  polynomial  becomes 

A  (s)  =  d(s)lm  ~  lmsn  +i t„Ims"~'  +  . .  .  +  d7lms  +  d,lm 

(14  c) 


t,  =  4/m 
Step  2 

Construct  a  matrix  polynomial  D(s)  using  A(s)  in  eqn.  14 
0(s>  =02„tlsJn  +  0JnsJ"-1  +  .  .  .  +  03$  +  D, 

=  FrAr(-s)A(s)f 

=  Fr{lms2n  +  (f„ 

+  )sJ"-2  +  . . .  +  £•,'£,  I F  ( 1 5) 

where  l)2ntl  -  FrF  =  R  is  a  weighting  matrix  to  be  determined. 

Step  i 

Solve  the  block  weighting  matrices  Q  and  R  from  eqns.  10  and  15 
in  terms  of  adjustable  parameters,  or 

Or*.,  =  F'F 

D2„  =  RA„  -a!,R  =  F’(F„  -£„')£ 


and  X  is  in  the  phase-variable  block  co-ordinate  and  consists  of  two 
block  vectors  ( X ,,  i  =  1,  2)  and  each  vector  Af,  consists  of  two  state 
variables  (xl  /t  /  =  1 ,  2,  /  =  1 ,  2).  The  state  equation  in  the  phase- 
variable  block  co-ordinate  is 


*1,1 

*1.2 

*2.1 

*2.2 

2622  1  1  -  11-8567  262.21 


0  005214  -  136-829 


101  -368 J  Ijcj  J 


y  1 

14-98 

95150 

0 

0 

yi 

.85-2 

124000 

0 

0 

=  012  +  A  !  RA;  Qu-A  fRA , 


=  F'lEiF.-  -EfE,  )F 

18-5671 

A,  = 

—  2622-1 

■  A2  = 

11-8567 

-  262-21 

=  On  +>1(764,  =  £'£/£,  F 

(16) 

.-  0-005214 

136-829 

0 

101-368 

Step  4 

Determine  the  required  block  weighting  matrices  Q  and  R  by  ad¬ 
justing  the  assigned  unknown  parameters  such  that  R  is  positive 
definite  and  Q  is  nonnegative  definite  symmetric. 

The  procedures  can  be  well  illustrated  by  the  following  gas-turbine 
example. 


An  illustrative  example 

Consider  the  following  linearised  two-shaft  gas-turbine 

J-IS 

—  I  268  -004528  1-498  951-5]  [a, 

1-002  -1-957  8-52  1240  z. 


14-98 

95150 

0 

C-  = 

0 

.  85-2 

124000.' 

0 

0 

(19  d) 

It  is  required  to  determine  two  optimal  block  controllers  for  the  gas- 
turbine  system  by  using 

(а)  assigned  weighting  matrices  Q  and  R  of  the  quadratic  performance 
index 

(б)  assigned  control  specifications. 

The  procedures  are  described  as  follows: 

(a)  Optimal-block-controller  design  via  assigned  weighting  matrices 
The  cost  function  of  the  state  equation  in  the  original  co-ordinate 
in  eqn.  17  is 


10  0 
I  0  100 


y=-|  | zTQz  +  uTRu]dt  (20) 

2  *0 

where  Q  =  lA  and  R  =l2  that  were  suggested  by  Tiwari  et  al. 15 
The  corresponding  cost  function  in  the  phase-variable  block  co¬ 
ordinate  is 


y  =  -  I  [XTQX  +  uTRu]dt 
2  Jo 


10  0  0 
0  10  0 


The  state  equation  in  eqn.  17  is  a  system  formulated  in  a  general 
co-ordinate.  To  apply  the  proposed  method,  a  block-linear- 
transformation  matrix  T  is  determined  from  the  newly  developed 
method  shown  in  Appendix  8.  The  block  linear  transformation  is 

z  =  nr  (is) 


where  R  =  l2  and 


=  TTrtT  — 


Oil  On 
O21  022 


7828-1776 

1 1941461-5  1 

185-671 

-  -521389 

11941461-5 

24436416630  , 

-26221 

13682-9 

185-671 

-26221  l 

100 

0 

—  -521389 

13682-9  , 

0 

10000 

From  eqn.  10  we  have 


'  14-98 

95150 

1 

1 

0 

o' 

[01  02- 

05]  =  [/2  -h  h] 

ot, 

Q?2 

aTr 

O2 

85-2 

124000 

1 

1 

0 

0 

02 

02*. 

0*2 

a[r 

18-5671 

-2622-1 

1 

10 

0 

0, 

02 

RA, 

ra2 

-0005214 

136-829 

( 

1 

0 

100 
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By  expanding  eqn.  22. 0(r)  in  eqn.  1 1  becomes 
0<s)  -  Dss4  +  D*s3  +  .  .  .  +  0, 

-  Rs4  +  liMj  ~A2'70sJ 


+  (RA,  +  A',R~Q22  A‘2RA:  Is1 
+  ((?i2  +  4 !  «42  -  y2i  -  4 Ira ,  is 
+  ((>ii  +  4  [ra  ,  )  -  02 


where 


Os  =  /2./>4 


0  -262-21 
262-21  0 


O, 


-2  03447  486-84 

486-84  -88755-96 


(23) 


02 


0  52440-924 

-  52440-924  0 


0.  = 


8172-92  11892776 

11892776  2-44433  x  lO10 


Performing  the  spectral  factorisation7  on  the  D(s)  gives 
A(s)  =  /2s2  +£,  s  +  E, 
where 


(24) 


E, 


17-2396  -261-906 
0-30451  576-845 


and  £| 


46-925  -  3929-922 

77-27215  156294-2 


From  eqn.  13  we  have  the  optimal  block  controllers  in  the  block 
co-ordinate  and  original  co-ordinate  as 


u  =  (. 4 ,  -£,  :  42  -£a)3r 
28-3576  —  1307-82 


5-3829 


0-3045 


1 77-2774  1  56157-4  0-304513  475  4  76 

=  M,  -£,  '•  A2  -£3]  r'z 


* 

(25a) 


-0-36296  0  279346  0-53829  0  003045 

I 

0-598572  0-795425  i  0  0304513  4-75476 


(25b) 

{b)  The  optimal-block-controller  design  via  assigned  control  specifi¬ 
cations 

The  design  goals  are  specified  as  follows. 

(i)  static  decoupling 

(ii)  final  values  of  the  unit-step  responses  are  unity 

(iii)  peak  time  tp  that  is  the  time  required  for  the  unit-step  res¬ 
ponse  to  reach  the  first  peak  of  the  overshoot  is  near  0-0 1  s 

(iv)  maximum  percentage  overshoot  is  less  than  10%. 


The  choices 
assigned  at 


eqri 


T 


-  225  tj  19843 


(21c) 


From  eqn.  26  l)(s)  can  be  determined  as 
/J(j)  =  Fl  L‘(— s)tk(s)F 

=  F'Fs 4  +  <2co2  -4|Jco ilF'Fs1  +  u*„FrF 
~  Ms4  +  <7 IA2  —  A  2R  )s3 

+  (RA  i  +  A[R-Q22  -  A[RA2)s 2 

+  ((>i2  AA'[RA2  -  Q‘,2 -AiRA,)s 

+  (Q„  +A]'RA,)  (28) 


For  simplicity,  let  (Jl2  -  Q2,  =  02  Fquating  the  matrix  coefficients 
of  the  same  power  of  eqn  28.  we  obtain  the  following  matrix 
equations: 

(а)  R  =-  FrF  (29 a) 

(б)  RA2  -  A[R  02  (296) 

(f)  RA,  +A{R  ~Q22  A(RA2  (2cj2, -4pul)FTF 

(29r ) 

(d)  A ,  RA  2  —  A  2  RA ,  =  ()2  (2 9d) 

(e)  <?..  +  A[RA,  =  w*FrF  (29 e) 


R  is  an  m  x  m  symmetric  and  positive-definite  matrix  which  has 
m(m  +  l)/2  unknown  elements  to  be  determined.  The  left-hand- 
side  matrices  in  eqns.  296  and  29 d  are  skew-symmetric  matrices. 
Expanding  the  matrix  equations  in  eqns.  296  and  29 d  results  in 
ni(m  —  1 )  simultaneous  equations  with  m(m  +  l)/2  unknown  vari¬ 
ables  in  R.  In  general,  there  are  an  infinite  number  of  solutions. 
However,  if  k  independent  simultaneous  equations  exist,  and 
k  <  m(m  +  l)/2,  then  we  can  assume  |m(m  +  1  )/2  —  *]  c  instants 
to  solve  k  unknown  variables  in  R.  The  choice  of  the  assigned  con¬ 
stants  in  R  is  a  design  freedom  and  a  certain  amount  of  experience 
is  helpful.  In  this  example,  we  assume  R„ ,  which  is  the  first  leading 
diagonal  element,  is  unity.  Thus  we  can  solve  for  R  and  F  in  eqn.  29 a 
as 


I 

2-92934 


2-92934 

51058-01562 


=  FrF 


where 


0-999916 

0-0129466 


5-737  x  10"s 
225-96 


(30) 


Note  that  R  is  a  positive-definite  matrix.  From  eqns.  30,  29 c  and 
29c  we  can  solve  for  Qn  and  Q22  as 


w*  -  345-55808 
2-929341oo*  +  77629-05 


2-929341ui*  +  77629-05 
51058-01  56  —  7-6070451  X  10*  J 

<31fl> 


To  reach  the  first  design  goal,  the  characteristic  matrix  polynomial  is 
defined  as 

A(s)  =  l2s2  +  E2s  a  E,  (26) 


Q„ 


4f‘wi  -2wA  -  140-5813 
2-929341(4{’u4  -  2w;i  -  2844-916 


where 


E> 


2*<o„  0 

0  2?a>„ 


and  E, 


cuS  0 

0  <o2 


2-929341  (4E’w^  -  2wJ)  -  2844-916 
51058-01 56(4f  2w)j  -  2a4)  -  52456131 7-4 


(316) 


£  (damping  ratio)  and  ui„  (undamped  natural  frequency)  are  unknown 
parameters  to  be  determined.  To  satisfy  the  third  design  goal  we  can 
estimate  w„  from  the  following  rule  of  thumb  in  designs14  as 


to*  —  —  —  —  300  rad/s  (27a) 

tp  0-01 

Also,  from  another  rule  of  thumb,14  we  can  estimate  $  to  meet  the 
fourth  design  goal  as 


f  a  -- 


In  Af„ 


In  0-1 
’  3-14 


—  0-75 


(276) 


Substituting  u>„  =  300  and  |  =  0-75  into  eqn.  31  yields  positive- 
definite  matrices  Q„  and  Q22.  Thus  the  optimal  block  controllers 
can  be  easily  found  in  the  block  co-ordinate  and  in  the  original 
co-ordinate  as 

u  =  [A,  —Et  •:A2-E1)X 


899814329  2622  1  1 

0-005214  89863-1 7  J 


438-1433  262-21  ' 
0  348-632. 


X  (32 a) 
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Ml  -e, ;  A3  -e2]t-'z 
—  1767-7  1357  .'7  1  43-8143  2-6221 

1-2182  -0-21391  ,  0  3-48632 


The  block-weighting  matrix  Q  in  the  block  co-ordinate  and  the 
weighting  matrix  Q in  the  original  co-ordinate  are 


8-1  x  10’  2-37277  x  1010 

2-37277  x  1010  4-13569  x  1014 

0  0 

0  0 


To  achieve  the  first  and  second  design  goals  we  sdd  a  forward-gain 
matrix  H  as  shown  in  Fig.  1 .  The  H  can  be  solved  from  the  block  C, 
ineqn.  I9ri  or 


14-98 

95150 

f-  198-42349 

152-258 

.85-2 

124000 

0-136336 

-0023971 

0  0 

0  0 

1 1703-42  31850-2 
31850-2  80169820  J' 

(33a) 


Thus  the  design  system  is 

•V,(l)]  =  1 
y2(s)  s2  +  2(co„s  +  cj2  [o 


o]  r  /?,< 

lj  [*,(. 


90000 

s2  +  450s +90000 


i  oj  r  *,( 

o  iJUj( 


3253160 

-2454610  1 

47-2383 

-  2-91217 

2454610 

1878440 

-33-808 

-  5-9383 

47-2383 

-  33-808  1 

117-0342 

31-850 

2-91217 

-  5-9383  | 

31-8502 

8016-982 

It  is  noticed  that2  any  arbitrarily  prescribed  closed-loop  poles  or 
control  specifications  may  not  result  in  a  positive-definite  matrix  R 
and  nonnegative-definite  matrix  Q.  The  constraints  suggested  by 
Anderson2  should  be  satisfied.  In  addition,  some  realistic  constraints 
to  the  amplitudes  of  the  control  signals,  for  example  the  limitations 
of  the  actuator  amplitude  and  rate  change  of  amplitude,  should  be 
also  examined. 


For  this  real  nontrivial  system  the  designed  system  is  not  only  static 
decoupling  but  also  complete  noninteracting,  and  the  final  values  of 
the  unit-step  responses  are  unity.  The  peak  time  is  0-014  s  and  the 
maximum  percentage  overshoot  is  1%.  The  simulation  curves  for  unit- 
step  input  are  shown  in  Figs.  2  and  3.  Comparing  the  design  results 
of  the  proposed  method  with  those  of  McMorran14  and  Tiwari  et 
al.,'s  the  present  result  gives  less  overshoot  and  less  oscillatory  res¬ 
ponses. 


Fig.  1 

Structure  of  designed  system 


Fig.  2 

Responses  of  various  designed  systems  to  a  unit  step  in  r , 

r,  =  I 

r2  =  0 

- McMorran’s  method: 

- proposed  method:  C_~  0*75;  u>n  =  300 

-  proposed  method:  Q  =  /4:/?  =  /, 

Tiwari’s  method:  Q  ~  /4;  R  =  /, 


$ 


b 


Fig.  3 

Responses  of  various  designed  systems  to  a  unit  step  in  r2 
r,  =  o 

r2  =  1 

- McMorran’s  method 

- proposed  method:  $  —0-75;  urn  —  300 

- proposed  method:  ~  /4  ;  R  =  I2 

Tiwari’s  method:  Q  =  /4 ;  R  = 


5  Conclusion 

A  new  method,  based  on  a  state  equation  in  the  phase- 
variable  block  co-ordinate,  has  been  presented  to  determine  the 
optimal  block  controllers  for  a  class  of  multivariable  systems.  The 
reverse  process  of  obtaining  the  optimal  block  controllers  has  been 
used  to  determine  the  weighting  matrices  of  the  quadratic  perfor¬ 
mance  index. 

When  a  multivariable  dynamic  system  is  formulated  in  a  matrix 
differential  equation,  the  proposed  method  is  more  suitable  for  the 
determination  of  the  optimal  controllers  than  the  conventional 
approach.  Also,  it  is  simpler  to  determine  the  weighting  matrices  than 
the  conventional  approaches.  However,  the  proposed  method  is  limited 
to  a  class  of  multivariable  systems  whose  state  equations  can  be 
formulated  into  matrix  differential  equations  or  the  state  equations 
in  the  block  companion  form. 
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matrices,  respectively.  matrix  transfer  function  or 

eqn.  38  can  be  directly  formulated  as 

YU)  =  [TV,  +jV2j  +  ...  +  N„ik- '|[0,  +D2s+  ... 

+  V'1  +/mJ*|-|fT(s) 

=  Ms)D-'(s)VU)  =  Tr(s)U(s)  (39) 

where  VU)  and  YU)  are  Laplace  Iransformsof  u(l)  and >>(/).  and  TrU) 
is  a  malrix  transfer  function. 

The  objective  is  to  derive  the  linear-transformation  malrix  7j  in 
eqn.  37.  Because  T ,  transforms  a  state  equation  in  eqn.  36  to  a  block 
companion  form  in  eqn.  38,  7",  is  called  as  a  block  linear  transform¬ 
ation.  We  further  assume  that  the  matrix  B0  in  eqn.  36  can  be 

_ 1  *L.-  r - f^lll  ...t _ n  rr  i.n.  X  m  ..  _ _ 


partitioned  into  the  form  of 


where  B2 1  £  R" 


singular  matrix.  This  can  be  accomplished  by  rearranging  the  sequence 
of  the  elements  in  the  state  vector  jr0 (/)  in  eqn.  36.  By  applying  the 
first  linear  transformation 

x„(t)  =  *,*,(/)  (40) 
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we  have 
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8  Appendix 

Block  linear  transformation 

Consider  a  class  of  completely  controllable,  linear,  time- 
mvariant.  multi-input,  multi-output  system 

i0(r)  =  A0x0(t)  +  B0u(r)  (36  a) 

yU)  =  c0x0(t )  (36b) 

where  A„£«"x",  Ba  £«"Xm,  C„  e  R1*"’  x„(t)  e  R"  x  1 , 
y(f)  E  R,x  1 ,  u(l)  £  Rm  x  1 .  Assume  that  l,m  <n  and  n/m  =  k  (an 
integer)  and  define  r  -  n  —  m.  By  a  linear  transformation 

x0U)  =  r,2,(r)  (37) 

We  wish  to  construct  a  state  equation  in  the  controllable  block 
companion  form 

z,(r)  =  A,z,(r) +«,u(r)  (38a) 

y(t)  =  C,z,(i)  (38b) 
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To  obtain  the  required  state  equation  in  eqn.  38.  we  perform  the 
second  linear  transformation 


(36a) 

x.  (/)  —  K2z. 

where 
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T  designates  the  transpose  of  the  matrix.  The  unknown  matrices 
Ql  £  Rr  x  r  (with  r  column  vectors  qs)  and  QT  ERrXm  (with  m 
column  vectors  ?,-)  can  be  evaluated  as  follows. 

From  eqn.  42a,  4 la  and  38c  wc  have  the  matrix  equation 

Kj'A,  =  A0K2]  (43a) 
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A, i  E  RrX  r,A  l2  E  Rr  *  m  ,A2.  E  Rm  x  r.  and  A22  ERmX  m . 

The  constant  matrices  0,  E  Rm  x  m  and  N,  €  R1  x  m  are  called  block 
elements  and  the  matrix  lm  =  lm  x  m  £  Bm  x  m  is  an  identity  matrix. 
The  matrices  Om  =  Om  x  m  €  Rm  x  m  and  Or  x  m  €  R’ x  m  are  null 


Expanding  eqn.  43b  yields 

0.4ii  =  AuQi  +-4|202 

0.-4. 2  =  -4  12  (43c) 

and 

02-4.1  +  A2i  —  -42,Q,  +  -422Q2 

02-4,2  +  A22  =  A  22  (43 <1) 

Performing  a  transpose  operation  on  eqn.  43c  and  substituting  eqns. 
38c  and  42c  into  it,  we  have  the  following  recursive  formulas: 

-4  =  qm,t  for  z  =  1,2 . r  (44 a) 

-4.2?t  —  Om  x  i  for  i  =  1, 2, . . .  ,r  -m  (44 b) 
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and 


and 


=  r1  for  /  =  1.2, (44c) 

where  e*  is  the  m  x  1  unit  column  vector  whose  ith  element  is  unity, 
and  all  other  elements  are  zeros.  Eqn.  44  can  be  further  simplified  as 
follows: 

(i)  If  k  =  2  then 

0,  =  Mil)"1*'  for  i  =  1,2,  ...,m  (45  a) 
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Applying  the  recursive  algorithm  in  eqn.  45 a,  we  have 
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When  the  square  matrices  in  eqns.  45a  and  45c  are  not  singular,  the 
in  eqn.  42  can  be  obtained.  Note  that  the  determination  of  qt  in 
eqn.  45  only  involves  one  inversion  of  a  matrix.  Thus  the  transform¬ 
ation  matrix  T ,  in  eqn.  37.  which  links  the  co-ordinates  x0(l)  in 
eqn.  36  and  the  required  co-ordinates  z,  (r)  in  eqn.  38,  is 

x0(r)  =  *.*21,(0  =  T,z,(r)  (46) 

It  is  believed  that  the  block  linear  transformation  T,  is  new. 

An  illustrative  example 

Consider  the  dynamic  equation  of  an  actual  gas-turbine  system13 
which  is  completely  controllable  and  observable. 

x0(t)  =  A0xa(t)  +  B„u(r) 

■HO  =  C0jr0(r)  (47) 

where 
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The  transformation  matrix  K2  in  eqn.  426  is 
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The  block  linear  transformation  T,  in  eqn.  46  is 
x0(t)  =  *,*22,(0  =  r,z,(  0 


h  =  4,  l  -  m-2,  r  =  n—  m  =  2,  and  k  =  n/m  =  2.  The  block  com¬ 
panion  form  in  eqn.  38,  the  corresponding  matrix  transfer  function, 
of  this  system  are  required. 

Applying  the  linear  transformation  in  eqn.  40  yields  the  state 
equation  in  eqn.  41 

i,(r)  =  A,*,  (r)  +  fl,u(r) 

y-(0=C,x,(r)  (48) 
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The  required  block  companion  form  in  eqn.  38  is 

r,(0  =  A,2,(r)  +  «,it(0 
.HO  =  C,z,(0 

where 
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The  corresponding  matrix  transfer  function  in  eqn.  39  is 
K(s)  =  [Ni  +  *2*1  (/),  +D2s+I2s2]-'U(s ) 
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Determination  of  Equivalent  Dominant  Poles  and 
Zeros  Using  Industrial  Specifications 

LEANG-SAN  SH1EH,  member,  ieee,  YING-JYI  PAUL  WEI,  member,  ieee.  HSI-ZEN  CHOW,  and  ROBERT  E.  YATES 


Aftwnsrf-A  graphical  method  and  an  analytical  method  are  presented 
to  determine  the  equivalent  dominant  poles  and  zero*  ot  a  system  using 
amigned  industrial  specifications.  A  second- order  tnnsfei  function 
with  two  poles  and  one  finite  zero  is  used  to  investigate  the  relation- 
riiipt  between  industrial  specifications  and  the  two  poles  and  one  finite 
zero.  Alto,  it  is  used  to  verify  the  tule  of  the  thumb  obtained  from 
Axelby’s  empirical  results.  A  frequency  response  data  matching  method 
is  proposed  for  fitting  a  low-order  transfer  function  using  the  assigned 
industrial  specifications  that  are  obtained  from  a  given  high-order 
transfer  function.  Thus  the  equivalent  dominant  poles  and  zeros  of 
sh%h  -order  system  can  be  determined  from  the  identified  low-order 
model. 


I.  Introduction 

IN  the  filter  and  compensator  designs  it  is  necessary  and  use¬ 
ful  to  have  a  rapid  method  or  a  simple  graphical  method  to 
determine  the  poles  and  zeros  that  dominate  the  characteris¬ 
tics  of  the  transient  response.  These  poles  and  zeros  are  called 
the  dominant  poles  and  zeros  that  can  be  used  to  estimate  the 
dynamic  behavior  of  the  system  response.  In  the  literature, 
the  definitions  of  the  dominant  poles  and  zeros  are  ambiguous. 
For  example,  the  dominant  poles  are  commonly  defined  as  the 
poles  which  are  located  near  the  imaginary  axis  (the  /w  axis) 
or  the  poles  which  have  the  smallest  absolute  value  when  no 
significant  zeros  appear.  Sometimes  a  pole  Pl  is  defined  as  the 
dominant  pole  [1]  if  l/><l>6l/*|l  where  Pt  are  other  system 
poles.  The  roles  of  dominant  zeros  that  are  often  neglected 
in  the  literature  become  significant  if  the  precise  dynamic 
characteristics  of  a  system  in  the  transient  state  are  required. 
The  zeros  not  only  contribute  to  the  initial  conditions  of  the 
transient  response  but  also  increase  the  bandwidth  in  the  fre¬ 
quency  domain;  therefore,  the  roles  of  the  zeros  are  as  im¬ 
portant  as  those  of  the  poles. 

As  the  technologies  are  progressing,  the  accurate  description 
of  many  physical  systems  results  in  a  high-order  transfer  func¬ 
tion  that  consists  of  many  clustery  poles  and  zeros  in  the  s 
plane.  The  poles  near  the  /go  axis  may  not  be  dominant  poles 
because  the  dominant  effects  on  the  transient  response  be¬ 
havior  of  the  poles  are  cancelled  by  the  nearby  zeros,  and  the 
system  response  may  be  characterized  by  the  collective  efforts 

Manuscript  received  July  13,  1978;  revised  January  25,  1979.  This 
work  was  supported  in  part  by  U.S.  Army  Missile  Command,  Redstone 
Arsenal,  AL.  DAAK  00-79-C-0061,  and  U.S.  Army  Research  Office 
DAAG29-77-G-0I43. 

L.  S.  Shieh,  Y.  J.  Wei,  and  H.  Z.  Chow  are  with  the  Department  of 
Electrical  Engineering,  University  of  Houston,  Houston,  TX  77004. 

R.  E.  Yates  is  with  the  Guidance  and  Control  Directorate,  U.S.  Army 
Missile  Research  and  Development  Command,  Redstone  Arsenal,  AL 
35809. 


of  a  group  of  clustery  poles  and  zeros.  This  implies  that  the 
poles  and  zeros  which  are  not  near  the  /«  axis  may  dominate 
the  characteristics  of  the  system  response.  Therefore,  the 
equivalent  dominant  poles  and  zeros,  rather  than  the  dominant 
poles  and  zeros  obtained  from  the  geometric  locations  in  the 
s  plane,  become  significant  in  the  analysis  and  synthesis  of  a 
high-order  system.  Furthermore,  the  design  goals  and  the 
nature  of  a  high-order  system  are  often  characterized  by  a  set 
of  control  specifications  [2]  (called  the  industrial  specifica¬ 
tions)  that  are  commonly  classified  as  1)  the  time-domain  spec¬ 
ifications,  for  example,  the  rise  time  and  the  overshoot,  2)  the 
frequency-domain  specifications,  for  example,  the  bandwidth 
and  the  phase  margin,  3)  the  complex-domain  specifications, 
for  example,  the  damping  ratio  and  the  undamped  natural 
angular  frequency  or  the  equivalent  poles  and  zeros  in  the  s 
plane.  If  the  relationships  among  the  time-domain,  frequency- 
domain  specifications,  and  the  equivalent  poles  and  zeros  (the 
complex-domain  specifications)  can  be  simply  determined 
from  a  simple  equation  or  a  working  graph,  then  the  selected 
poles  and  zeros  in  the  design  of  filters  and  compensators  be¬ 
come  meaningful,  and  the  design  processes  can  be  greatly 
simplified. 

In  this  paper,  a  graphical  method  and  an  analytical  method 
are  proposed  to  determine  the  equivalent  dominant  poles  and 
zeros  using  assigned  industrial  specifications.  First,  relation¬ 
ships  among  various  industrial  specifications  will  be  studied. 
A  second-order  transfer  function  having  two  poles  and  one 
finite  zero  is  used  as  a  basis  for  the  investigation.  Several 
working  graphs  and  mathematical  expressions  are  developed 
for  the  determination  of  the  two  dominant  poles  and  one 
dominant  zero  using  the  assigned  industrial  specifications.  Then 
the  equivalent  dominant  poles  and  zeros  of  a  high-order  sys¬ 
tem  are  determined  by  a  new  dominant  frequency-response 
data  matching  method.  The  equivalent  dominant  poles  and 
zeros  thus  obtained  satisfy  the  exact  assigned  industrial 
specifications. 

II.  The  Relationships  Among  Various  Industrial 
Specifications 

In  control  system  design,  the  design  goals  are  usually  ex¬ 
pressed  in  terms  of  a  set  of  industrial  specifications.  The  place¬ 
ment  of  poles  and  zeros  based  upon  the  assigned  specifications 
needs  certain  experiences.  If  the  relationships  among  various 
industrial  specifications  can  be  determined,  then  nonconflict¬ 
ing  industrial  specifications  can  be  assigned  as  design  goals,  and 
the  meaningful  dominant  poles  and  zeros  can  be  selected  for 
filter  and  compensator  designs.  Thus  an  effective  design 
method  may  be  developed. 
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An  empirical  study  on  the  relationships  among  various  in¬ 
dustrial  specifications  has  been  conducted  by  Axelby  [3) . 
The  empirical  rules  or  the  rule  of  the  thumb,  which  link  the 
specifications  in  both  time  and  frequency  domains,  are  listed 
as  follows: 


M,^Mp 
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sin  <p„ 


(la) 


where  M ,  is  the  maximum  value  of  unit-step  response,  Mp  is 
the  maximum  value  of  the  closed-loop  frequency  response, 
and  0m  is  the  phase  margin; 
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where  t„  is  the  peak  value  time  or  the  time  when  M,  occurs; 
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where  r„  is  the  time  when  the  maximum  error  of  the  ramp 
function  with  respect  to  its  input  occurs; 
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where  tc  is  the  time  when  Mt  occurs. 

Other  rules  of  the  thumb  according  to  Truxal  (4]  are  listed 
as  follows: 


tr co6  =  0.6rr  to  0.9rr 


(lh) 


where  tr  is  the  rise  time  or  the  time  required  for  the  response 
to  go  from  10  to  90  percent  of  its  final  value  and  ojb  is  the 
bandwidth  in  rad/s; 
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where  td  is  the  delay  time  or  the  time  required  to  reach  50 
percent  of  its  final  value  and  Kv  is  the  velocity  error  constant. 

Some  other  analytical  results  that  represent  the  relationships 
between  the  time-domain  specifications  (but  not  the  frequency- 
domain  specifications)  and  the  complex-domain  specifications 
have  been  developed  and  can  be  found  in  standard  textbooks 
[5] ,  [6) .  The  most  commonly  used  function  for  investigating 
the  relationships  is 


HS)  UJfi 

R(s)  s2  +  2l-tJ„s  +  <Vi 
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where  Y(s)  and  R(s)  are  the  output  and  input  functions, 
respectively,  and  £  is  the  damping  ratio  and  w„  is  the  un¬ 
damped  natural  angular  frequency.  From  (2)  we  observe  that 
the  zero  of  the  system  is  located  at  infinity,  and  is  not  a 
significant  zero.  Since  the  time-domain  specifications  are  of¬ 
ten  used  to  define  the  characteristics  of  the  transient  behavior, 
the  roles  of  zeros  become  significant.  Therefore,  a  better 
model  than  that  of  (2)  should  be  used  to  study  the  relation¬ 
ships  among  the  industrial  specifications.  The  transfer  func¬ 
tion  of  a  unit-feedback  system  that  has  two  poles  and  one 
finite  zero  is  used  as  a  basis  for  the  investigation.  The  pro¬ 
posed  closed-loop  transfer  function  is  then, 

,2 


where  Me  is  the  maximum  value  of  the  error  of  the  unit-ramp 
function  and  uic  is  the  gain  crossover  frequency; 

(1c) 

where  cup  is  the  peak  value  frequency  or  the  frequency  when 
Mp  occurs; 

M,^ioc  (Id) 

where  M,  is  the  maximum  value  of  the  unit-impulse  response; 
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where  s*  is  a  normalized  complex  variable,  a,  and  £>,  are  con¬ 
stants.  and  A(s)  and  B(s)  are  two  polynomials.  The  normalized 
poles  and  the  original  poles  are  at 
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and  the  normalized  zero  and  the  original  zero  are  at 
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The  open-loop  transfer  function  (/( s)  of  the  system  in  (3)  is 

where  Ku  =  co„/(2£  r)  is  the  velocity  error  constant  if  t  <  2% 
a  =  (2^  r)  to,,  and  b  =  t j„/t. 

Comparing  (2)  and  (3)  we  observe  that  a  finite  zero  has  been 
inserted  in  (3).  The  zero  contributes  the  initial  condition  at 
the  transient  state,  and  it  reduces  the  velocity  error  at  the 
steady  state.  Also  it  provides  an  additional  bandwidth  in  the 
frequency  domain,  which  increases  the  phase  margin  and  im¬ 
proves  the  stability  of  a  system. 

The  derivations  of  the  relationships  among  the  industrial 
specifications  arc  shown  as  the  following  seven  relationships. 

1)  The  Relationships  Among  M,.  tp,  £.  u>„.anj  t:  The  unit- 
step  response  of  the  system  in  (3)  gives 


TCJ»5  +  C0« 

T(s)  =  — — — — - 

s(sJ  +  2£t o„s  +  )’ 

The  inverse  Laplace  transform  of  )'(s)  results  in 
Y(t)  =  1  I 


(6a) 


cosa>„Vl  %2r 


f  r 

Vi  t: 


sin  u>„  vry 


(6b) 
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Differentiating  y{t)  with  respect  to  t  and  setting  the  result 
equal  to  zero  yields 


Substituting  (6c)  into  (6b)  and  simplifying  it  gives  the  maxi¬ 
mum  value  of  the  unit-step  response 

Af,*l+  -  2 +  1)1/2.  (6d) 


2)  The  Relationships  Among  Mp,  up,  cj„,  and  r:  Apply¬ 
ing  Higgins  and  Siegel’s  complex  variable  differentiation 
method  [7] ,  we  can  solve  the  peak  value  frequency  up  from 
the  following  equation: 


Re 


1  d£js) 
B(s)  ds 


1  <^(*)1 1 
A(s)  ds  J|,./u>p~ 


Thus  we  have 

cjp  =  u„  Vl  -  2(J 

MpM/(2iVrT)J’ 

and 


ifr  =  0 


(7a) 


(7b) 


<op  =  ^  (- 1  +  V(rJ  +  1)J  -  4tj{j  ] 1/1 

,2 


Mp  Mr1  +  l)1  -  dtV  -  (r1  +  1) 


>  ,  if  r  ^  0. 


+  2{V]*1/S- 

5)  77ie  Relationships  Among  <t>m,  cjc, 
the  definitions  of  0m  and  cjc, 

+  180° 

and 


J  (7c) 

{,  oj„,  and  r:  Using 

(8a) 


Ws-A*  *  1  (8b) 

we  have 


1  tan* 


1  -  (2 t-r)r 


and 

wc  3  «,  [2{r  -  2|a  +  V(2{*  -  2«r)*  +  1]  . 


(8c) 

(8d) 


4)  The  Relationships  Among  tv>  Me,  {,  u>„,  and  t:  The  er¬ 
ror  signal  e(t),  which  is  the  difference  between  the  ramp  in¬ 
put  r(t)  and  the  time  response  y(r)  of  the  same  input  to  the 
system  in  (3),  is 


e(t)  m  — — -  -  e**u*,[f?  cos  u„  Vl  -  $af 

wb  Au„ 

-  Csin  wBVl  *  f*fj  (9a) 

where 

A-(l-*J),*-(2{-rXl -tJ), 

0(1  -  2{a  +  r{)Vl  -  *a. 


Differentiating  e(t)  with  respect  to  t  and  setting  the  result 
equal  to  zero  we  have 


Substituting  the  tu  into  (9a)  and  simplifying  it  we  have 


Me  =  [2*  -  r  +  V(1  +  t5  -  2r{)e*  ^  ]  /«„ .  (9c) 

5)  The  Relationships  Among  tc,  M„  (,  un,  and  t:  Dif¬ 
ferentiating  the  unit-impulse  response  y(f)  of  the  system  in 
(3),  y(t),  and  setting  the  result  equal  to  zero,  we  have  the  time 
rc  at  which  the  maximum  value  occurs,  or 


■  tan’ 


r(l-2tr)Vn?1 
L  {-2t^+t  I 


(10a) 


Substituting  tc  intoy(r)  yields  the  maximum  value  of  the  unit- 
impulse  response  Mt,  or 

M,  =  w„e*^'cVr2  -  2{t  +  1.  (10b) 

6)  The  Relationships  Among  Kv,  co„,and  r :  The  velocity 
error  constant  Kv  can  be  derived  from  the  basic  definition  as 

JC.-Bma-GW-rJ5*-,  ifr<2{.  (11) 

i«o  2{-r 


7)  The  Relationships  Among  tob,  f,  and  t:  The  defini¬ 
tion  of  the  bandwidth  of  a  system  is 


(Tm-M,  -TTf 


'sp¬ 


in) 
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The  analytical  expression  is 

w»=wn[(l+f:  2£2)  +  V(  I  +  t2  27*7+71 11 2 .  (13) 

Most  important  time-domain  and  frequency-domain  specifi¬ 
cations  have  been  analytically  expressed  in  terms  of  j|,  u>„.  and 


lig.  4  Relationships  among  1/sin  e>,„.  £.  uj„.  and  r  shown  in  (8c). 


I  ig  5  Relationships  among  wt..  £.  w„.  and  t  shown  in  (8d). 

r  which  ate  the  specifications  in  the  complex  domain.  These 
expressions  arc  normalized  and  graphically  shown  in  Figs.  1- 
11.  If  an  industrial  specification  is  assigned,  the  corresponding 
£  and  r  or  the  equivalent  poles  and  zero  in  (4)  can  be  deter¬ 
mined  from  the  plotted  curves.  Also  the  curves  in  Figs.  12- 
15  can  be  used  to  verify  the  rules  of  the  thumb  proposed  by 
Axelby  [3],  It  is  observed  that  the  accuracy  of  the  rules  de- 
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Fig.  6.  Relationships  among  wp,  f,  u and  r  shown  in  (7). 


0.  0.25  0.5  0.75  '1.0 

Fig.  7.  Relationships  among  t„,  t,  wn,  and  r  shown  in  (9b). 


pends  upon  the  range  of  the  damping  ratio  and  the  zero  loca¬ 
tion.  Furthermore,  from  the  developed  working  graphs,  a  set 
of  meaningful  and  nonconflicting  specifications  can  be  as¬ 
signed  for  the  design  goals  of  a  control  system. 

III.  Determination  of  Equivalent  Dominant  Poles 
and  Zeros  from  a  High-Order  Model 

In  the  design  of  high  performance  control  systems,  quite 
often  several  specifications  are  assigned  as  design  goals,  and  the 
corresponding  dominant  poles  and  zeros  are  required.  This  is 
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0.  0.25  0.5  0.75  '.-1.0 


Fig.  8.  Relationships  among Me,  (,  u„,  and  r  shown  in  (9c). 


0.  0.25  0.5  0.75  C-1.0 


Fig.  9.  Relationships  among  Mt,  t,  a>„,  and  r  shown  in  (10b). 

a  problem  of  a  high- order  transfer  function  fitting  using  indus¬ 
trial  specifications.  Shieh  et  al.  [8J ,  [9]  have  developed  an 
original  synthesis  technique  to  fit  a  second-order  transfer  func¬ 
tion  based  on  three  industrial  specifications.  The  Newton- 
Raphson  multidimensional  method  [10]  was  applied  to  solve 
the  resulting  nonlinear  simultaneous  equations  that  can  be 
converted  to  a  single  variable  quadratic  equation.  However,  it 
is  well  known  that  the  Newton-Raphson  method  will  only  con- 


TreT 


Fig.  10.  Relationships  among  tc ,  f,  w„.  and  t  shown  in  (10a). 

verge  for  a  small  range  of  starting  values  or  the  initial  guesses. 
It  is  also  known  that  high- order  nonlinear  equations  have 
many  solutions  that  depend  heavily  on  the  initial  guess  used. 
For  general  nonlinear  equations  that  cannot  be  converted  to  a 
single  variable  equation,  the  Newton-Raphson  numerical 
method  may  not  converge  to  the  desired  solution  using  arbi¬ 
trary  initial  guesses.  In  this  paper,  the  original  synthesis  method 
[8] ,  [  9]  is  extended  for  modeling  a  high-order  transfer  func¬ 
tion  using  many  industrial  specifications;  and  an  analytical 
method  is  proposed  for  the  estimation  of  the  good  starting 
values.  Thus  the  desired  dominant  poles  and  zeros  can  be  de¬ 
termined  from  the  identified  transfer  function.  The  method 
can  be  well  illustrated  using  the  following  example. 

Suppose  that  the  poles  and  zeros  that  represent  the  follow¬ 
ing  given  industrial  specifications  are  required  to  be  determined. 

Type  “  1”  system  (14a) 

«c  the  gain  crossover  frequency  =  4.7  rad/s  (14b) 

0m  the  phase  margin  =  45.6°  (14c) 

Mp  the  maximum  value  of  the  closed-loop  frequency 

response  =1.5  (14d) 

c Op  the  peak  value  frequency  ■  3.5  rad/s  (14e) 

u>b  the  bandwidth  of  the  closed-loop  frequency 

response  =  6.5  rad/s.  (14f) 

The  assignments  of  the  specifications  in  (14)  closely  follow  the 
rules  shown  in  (1).  Therefore,  the  conflicted  assignments  can 
be  avoided.  The  first  two  are  the  open-loop  specifications, 
while  the  others  are  the  closed-loop  ones.  Three  equivalent 
poles  and  two  equivalent  zeros  that  represent  the  assigned 


Fig.  11.  Relationships  among  uib,  (,  u>„,  and  r  shown  in  (13). 


specifications  in  (14)  can  be  determined.  The  third-order 
model  is 

r(»)  _  n  ,  _  *(*  +  ZiXs  +  *s) 

R(s)  S  (sJ  +  2$o)ns  +  wj  X*  +  P) 

h,sJ  +  bjs  +  b3  , 

=  - 1 -  (15a) 

s3  +0|SJ  +a2s  +  a3 

where  K,  p,  £,  u>n,  zx,  and  or  the  corresponding  a,  and  b/ 
are  unknown  constants  to  be  determined.  Because  the  system 
is  a  type  “1”  system,  the  final  value  of  the  unit-step  response 
of  the  system  in  (15a)  is  unity  or 


r(0|r—  =  lim  s  /?(s)T(s) 

J-0 

s*o  ^sMs'+atS'+ajs  +  flj/  a3 


or  a)=b3. 

As  a  result,  (15a)  can  be  simplified  as 

W  ,  b,!3  +b,s+ a, 

R(s)  s3  +ats2  +a}s  +  a} 

The  open-loop  transfer  function  G(s)  is 

G(  )  m  b^  +  bis  +  at 

’  s[sJ  +  («,  -h,)s  +  (fl,  -  ft,)]' 
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o.  c.;'i  u.‘  0.75  -i.o 

Fig.  12.  Relationships  among  Mp,  M,,  and  1/sin  0m  shown  in  ( X). 


Fig.  13.  Relationships  among  Me,  tc,  and  l/wc  shown  in  (1). 


Following  the  definitions  shown  in  (14),  we  can  construct 
a  set  of  nonlinear  equations  /,(a i,  a2,  a3,  bt,  b2)  -  0  for  i  = 
1,2,-  -,5. 

The  definition  of  cot.  is 

\G(jojc)\=  1.  (16a) 

The  corresponding  nonlinear  equation  is 

/1(a1,a2,(2j,61,62)  =  (fl1  -  6,)2 u*  +  [to*  -  (a2  *2)wc)2 

-  (a;  -  6,co2)2  -  h2u>2  =  0.  (16b) 

The  definition  of  <t>„.  can  be  expressed  as 
<Pm  ~  180°  +/G(juc)  .  (17a) 

The  nonlinear  equation  is 

/2(a|,a2,a3,6|,62)  =  h2a)2(a1  bx) 

*  (a3  6,w?Xtj2  a2  +  b2) 
tan  <t>,„  [(a  3  hico2)(a,  ~  b ,  )ojc 
+  62wf(w2  -  a2  +  62)J  =  0.  (17b) 

The  definition  of  u>b  is  known  as 

\TXiuh)\  =  ^=.  (18a) 

The  corresponding  nonlinear  equation  is 
/3(a,,a2.a3,61,62)  =  (a3  -  b, to£)2  +  b\u\ 

'  jl(o3  + 

=  0.  (18b) 


The  definition  of  u>p  gives 
rfl7\/'w)l| 

-r-H  =  0  (19a) 

doi  Icj  =  t Op 

Following  Higgins  and  Siegel's  complex  variable  differential 
technique  [7] ,  we  have  the  following  nonlinear  equation: 

Mai, a2,a3,6i,f>2)=  [2a ia3Gjp  -  2a2 wj 

-  (a3  -  3w2X-  Wp  +  «2tOp)]  [(a, 

-  bx co2)2  +(h2ojp)J]  +  [  2a36,wp 
+  2 b]u>3p  +  b\ Wp]  ((a3  -  a,Wp)2 


+  (-co2 +a2Op)2]  =0. 


The  definition  of  Mp  is 


|7TM)I 


Ul»Wp  ~Mp. 


The  nonlinear  equation  is 


/5(a,,a2,a3,6,,62)  =  (a3  -  6,w£)2  +hfw2 
-Af2[(a3  -  a,w2)2 

+  (w2  a2wp)2l  =0.  (20b) 

Equations  (16)-(20)  are  a  set  of  high- order  nonlinear  simulta¬ 
neous  equations  which  are  very  difficult  to  solve.  The  Newton- 
Raphson  method,  which  is  available  in  most  digital  computers 
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[ID  .  is  used  to  solve  the  nonlinear  equations.  To  obtain  the 
desired  solution,  and  to  improve  the  speed  of  convergence  of 
the  numerical  method,  we  have  to  establish  a  set  of  good  start¬ 
ing  values.  From  the  developed  analytical  expressions  of  vari¬ 
ous  specifications  or  the  working  curves  in  this  paper,  we  can 
determine  the  corresponding  two  poles  and  one  zero  using 
Mp  =  1.5  and  cop  =  3.5.  From  the  rule  of  the  thumb  in  ( 1 )  we 
observe  that  the  Mp  and  up  have  indirectly  included  the  ap¬ 
proximated  respective  0„,  and  cu, .  The  procedures  are  shown 
in  the  following  steps. 

Step  l:  Determine  the  normalized  dominant  poles  or  the  f  in 
(4a)  using  the  curve  drawn  in  Fig.  2.  having  r  =  0.  From  the 
curve  (t  =  0)  we  read  the  damping  ratio  %  =  0.35.  The  normal¬ 
ized  dominant  poles  and  the  dominant  poles  with  u>p  =  uj„  = 
3.5  arc- 

si^  0.35  +/0.9368  s,  =  1.225  + /3.278ft 

s  J=  0.35  /0.9368  s2  =  1.225  /3.278b.  (21a) 

The  second-order  model  is 


77(s)  = 


s2  +  2?w„s  +  s2  r  2.45s  +  1  2.25 


Step  2:  Determine  a  dominant  zero  using  the  specification 
u> h  =  6.5  in  (14f).  The  modified  second-order  model  becomes 

bts  +  bts  +  1 2.25 

T  = - - - - - = - - - ni,.i 

1  jJ  +  2|oj„s  +  coj  s2  +  2.45s  +  12.25 

The  b  i  can  be  easily  determined  by  using  the  definition  of 
Wi,  in  (18a): 

b,  -  3.1781.  (2 Id) 


Thus  a  low-order  dominant  model  is  determined.  However,  a 
third-order  model  is  required.  An  extra  pole  and  a  nearby 
zero  are  inserted  into  the  second-order  model  in  (21c)  to  ob¬ 
tain  an  approximate  third-order  model,  or 

t-*,  ._  <(>is  +  cj;,KI.1s  +  1 0i;U>„) 

7  ns)-  ,  -  ■ ,  . 

(V  +  2£u ;,,s  +  cc,‘,)(s  +  10^u!„) 

_  3.495‘Mr  +  52.406725s  +  150.0625 

s3  +  14.7s2  +  42.2625s  +  150.0625  '  (“'e) 

Using  the  coefficients  in  (21c)  as  initial  guesses:  a*  =  14.7 .a*  - 
42.2625.  at  =  150.0625.  I>t  =  3.49591.  and  6?  =  52.406725. 
and  applying  the  Newton-Raphson  method  [11]  to  solve  the 
nonlinear  equations  in  (16)  through  (20)  yields  the  desired 
solutions:  a,  =  4.267162.  a2  =  20.58799. a3  =29.806197 .b,  = 
3.188355.  and  h2  =  1 5.561058.  at  10th  iteration  with  the  er¬ 
ror  tolerance  of  10"1’.  The  desired  transfer  function  is 

3. 1 88355.v:  +  15.561058s  +  29.806197 

7"3(s)  =  -j - , - .  (22) 

s3  +  4.267162  s2  +  20.58799s  +  29.806197 

The  dominant  poles  and  zeros,  which  represent  the  assigned  in¬ 
dustrial  specifications,  arc  determined  from  the  poles  P ,  and 
zeros  r,  in  (22); 

=  1.849412756 

P2  =  1.208824622  +/3.8282263I8 

=  1.208824622  /  3.828226318  (23a) 


=  4.880591402  +  /3 .684243  78 
r2  =  4.880591402  /3.68424378. 
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When  the  distribution  of  the  poles  and  zeros  of  a  high-order 
transfer  function  is  known  and  the  reduced-order  transfer 
function  that  consists  of  equivalent  dominant  poles  and  zeros 
is  required,  it  is  a  model  reduction  problem.  Recently,  various 
model  reduction  methods  [12] -[IS]  have  been  proposed  in 
the  frequency  domain.  However,  their  reduced  models  [12]- 
[15] ,  do  not  keep  the  assigned  industrial  specifications,  which 
are  obtained  from  the  original  system.  The  preservation  of  the 
exact  frequency-domain  specifications  is  essential  in  the  de¬ 
sign  of  filters  and  compensators  using  frequency-domain  meth¬ 
ods  [S],  [6],  such  as  the  Nyquist,  Bode,  and  Nichols  chart 
methods.  This  proposed  method  can  overcome  the  short¬ 
comings  of  the  existing  model  reduction  methods.  The 
frequency -response  data  at  u>p,  uc,  and  (the  phase 
crossover  frequency  of  the  open-loop  system  for  the  use  of  the 
gain  margin  [5],  [6]  are  considered  as  the  dominant  frequency- 
response  data.  If  some  of  these  data  are  assigned  to  determine 
the  corresponding  reduced-order  model,  the  equivalent  domi¬ 
nant  poles  and  zeros  can  be  determined  from  the  reduced- 
order  model  that  consists  of  the  exact  industrial  specifications 
assigned. 

IV.  Conclusion 

A  second-order  transfer  function  with  two  poles  and  one 
finite  zero  has  been  used  to  derive  the  analytical  and  graphical 
expressions  of  various  industrial  specifications.  For  a  few  as¬ 
signed  industrial  specifications,  the  corresponding  two  domi¬ 
nant  poles  and  one  dominant  zero  can  be  determined  from 
the  identified  transfer  function.  The  generalized  second- 
order  model  has  been  used  to  verify  the  rule  of  the  thumb 
proposed  by  Axelby.  It  has  been  observed  that  the  accuracy 
of  the  rule  of  the  thumb  depends  on  the  range  of  the  damping 
ratio  and  the  zero  location.  From  the  developed  graphical 
expressions,  a  set  of  meaningful  industrial  specifications  can 
be  chosen  and  assigned  as  the  design  goals  for  the  filter  and 
compensator  designs.  A  dominant  frequency-response  data 
matching  method  has  been  developed  to  construct  a  low- 
order  transfer  function  using  the  assigned  industrial  specifica¬ 
tions  that  are  obtained  from  a  given  high-order  system.  Thus 
the  equivalent  dominant  poles  and  zeros  of  a  high-order  sys¬ 
tem  can  be  determined  from  the  identified  low-order  transfer 
function  that  has  the  exact  industrial  specifications  assigned. 


More  over,  the  proposed  method  in  this  paper  has  been 
successfully  applied  to  redesign  the  compensators  of  a  stabi¬ 
lized  pitch  control  system  of  a  real  semiactive  terminal  homing 
missile  [16].  The  overall  system  characteristics  of  the  rede¬ 
signed  missile  [17]  match  those  of  the  lower  ordered  model 
obtained  from  assigned  industrial  specifications. 
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A  method  for  modelling  transfer  functions  using  dominant 
frequency-response  data  and  its  applications 

L.  S.  SHIEHf,  M.  DATTA-BARUAf,  and  R.  E.  YATES  J 


This  paper  presents  a  fundamental  method  for  modelling  transfer  functions  using  the 
basic  performance  specifications  and  frequency -response  dai  .  at  the  dominant 
frequencies.  A  set  of  non-linear  equations  is  constructed  from  the  definitions  of  the 
basic  performance  specifications,  the  dominant  frequency -respond  data  and  the 
unknown  coefficients  of  a  transfer  function.  A  Newton-  Raphson  multidimensional 
method  is  applied  to  solve  the  non-linear  equations.  F>  ir  methods  are  given  to 
construct  approximate  representations  of  the  desired  transfer  functions  for  the 
estimation  of  guod  starting  values  to  ensure  rapid  convergence  of  the  numerical 
method.  The  applications  of  the  proposed  method  are  :  (l)  developing  a  standard 
model  and/or  a  transfer  function  of  a  filter  or  a  compensator  using  the  specified 
dominant  frequency -response  data  ;  (2)  identifying  the  transfer  function  of  a  system 
from  available  experimental  frequency -response  data  ;  and  (3)  reducing  high-order 
transfer  functions  to  low-order  models  using  dominant  frequency -response  data. 


1.  Introduction 

The  nature  of  the  transient  response  of  a  system  is  often  characterized  by  a 
set  of  performance  specifications  in  the  time  domain  such  as  the  settling  time 
and  the  rising  time.  In  the  frequency  domain,  another  set  of  performance 
specifications  (Gibson  and  Rekasius  1961)  is  used  to  represent  the  charac¬ 
teristics  of  the  system  performance.  The  bandwidth  and  the  phase  margin 
are  typical  examples  of  the  frequency  domain  specifications.  In  designing 
compensators  and  filters,  and  in  predicting  the  nature  of  time  response  of  a 
system,  practicing  engineers  are  often  interested  in  the  dominant  poles.  These 
can  be  converted  to  a  damping  ratio  and  a  natural  angular  frequency  specified 
in  the  complex  plane.  These  specifications  are  often  called  the  complex-domain 
specifications.  The  engineer  is  also  interested  in  various  error  constants  (for 
example,  the  velocity-error  constant),  which  represent  the  characteristics  of 
system  performance  in  both  time  and  frequency  domains  (Truxal  1955).  The 
frequency-response  data  at  the  frequencies  of  the  frequency-domain  specifica¬ 
tion  are  considered  as  the  dominant  frequency-response  data  in  this  paper 
because  these  data  characterize  the  nature  of  the  system  responses.  For 
example,  the  phase  margin  (<£m)  of  a  system  at  the  gain -crossover  frequency 
(coc)  is  often  used  as  a  measure  of  additional  phase  lag  required  to  bring  the 
system  to  the  verge  of  instability.  Also,  if  the  phase  angle  of  the  open-loop 
system  at  the  wc  is  near  —  180°,  then  the  response  of  the  closed-loop  system  will 
be  oscillatory. 
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In  the  design  of  a  control  system  in  the  frequency  domain,  the  specifications 
discussed  above  or  the  dominant  frequency -response  data  are  usually  considered 
as  design  goals.  Various  frequency -domain  or  complex-domain  approaches 
(Nyquist  1932,  Evans  1953,  Bode  1954,  Thaler  1973)  have  been  developed  and 
widely  applied  in  industry  for  compensator  designs  to  achieve  desired  perfor¬ 
mance.  The  most  popular  design  methods  are  those  based  on  the  Nyquist 
(1932)  plot,  the  Bode  (1954)  design,  and  the  root-locus  method  (Evans  1953, 
Thaler  1973).  To  improve  the  efficiency  of  the  design  Methods,  it  is  advan¬ 
tageous  to  have  the  design  goals  expressed  as  mathematical  functions  or 
transfer  functions  (defined  as  the  standard  models).  Once  standard  models 
have  been  ascertained,  the  corresponding  time-domain  specifications  and 
temporal  responses  can  be  determined  from  digital  or  analogue  simulations  of 
the  standard  models.  Also,  the  frequency-response  data  of  the  desired  com¬ 
pensator  can  be  determined  from  Nyquist  plots  or  Bode  diagrams  by  comparing 
the  frequency-response  curves  of  the  original  and  the  desired  response  models. 
The  required  filters  and  compensators  (Del  Toro  and  Parker  1900,  Thaler  1973) 
can  then  be  easily  determined. 

Empirical  rules  or  rules  of  the  thumb  that  link  the  specifications  in  the 
time,  frequency,  and  complex  domains  have  been  developed  by  Truxal  (1955), 
Del  Toro  and  Parker  (1960),  Axelby  (1960),  and  Seshadri  (1969)  et  al.  From 
these  results,  it  is  observed  that  most  time-domain  specifications  and  complex- 
domain  specifications  can  be  approximately  converted  to  frequency-domain 
specifications.  Some  of  these  frequency-domain  specifications  are  phase 
margin  maximum  value  of  the  closed-loop  frequency  response  (d/p), 

gain-crossover  frequency  (wc.),  peak  value  frequency  (tup),  the  bandwidth  (cub), 
and  velocity -error  constant  (A'v).  Other  important  frequency-response  data 
are  :  (I)  the  real  part  of  the  open-loop  transfer  function  G(j to)  at  the  phase- 
crossover  frequency  (c o„)  which  has  been  used  to  define  the  gain  margin  (Gm) ; 
(2)  the  real  part  and  imaginary  part  of  the  closed-loop  function  (7’(.j))  and  the 
open-loop  function  G(s)  at  s  =  ja>  —jajQ  =  JO.  The  data  at  w  =  0  often  indicate 
the  final  value  and  the  type  of  the  system.  In  a  type  1  system,  /m[G'(J0)]  has 
an  infinite  value,  while  Re  [G(j0)]  has  a  finite  value  from  which  an  asymptotic 
line  (Del  Toro  and  Parker  1960)  can  be  drawn  in  a  Nyquist  plot ;  (3)  the  corner 
frequencies  in  the  Bode  plot  of  G{ju>)  in  the  regions  of  w  =  wcl  where 
20  log  |G(Jajcl)|  =  +  15  dB,  and  oj  =  tue2  where  20  log  |G(jcuc2)|  =  -  15  dB. 
Chen  (1957)  has  shown  empirically  that  the  open-loop  poles  and  zeros  of  a 
system  can  be  approximated  by  retaining  the  Bode  plot  in  the  regions  of  the 
±  15  dB  boundaries.  Some  dominant  frequency-response  data  are  indicated 
in  Fig.  1. 

Through  use  of  the  above  dominant  frequency -response  data,  a  basic  method 
is  proposed  in  this  paper  for  modelling  various  transfer  functions.  First,  a  set  of 
simultaneous  non-linear  algebraic  equations,  based  on  basic  definitions  of  the 
dominant  frequency-response  data  and  the  unknown  coefficients  of  a  desired 
transfer  function,  is  constructed.  Then  the  Newton-Raphson  method 
(Carnahan  et  al.  1969,  IBM  1977)  is  used  to  solve  the  non-linear  equations. 
However,  as  is  well  known,  the  Newton-Raphson  method  will  often  only 
converge  for  a  small  range  of  starting  values  ;  therefore,  four  methods  are 
developed  in  this  paper  for  estimating  good  starting  values  so  that  the  numerical 
method  (IBM  1977)  will  converge  rapidly  to  the  desired  solution. 
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Figure  1 .  Nyquist  plot  of  an  open-loop  system  0(s). 


The  applications  of  this  method  can  be  classified  as  follows. 

(1)  When  the  design  goals  are  predescribed  by  the  dominant  frequency  - 
response  data,  which  may  be  obtained  from  the  frequency-domain 
specifications  (Gibson  and  Rekasius  1961)  or  equivalent  ones  (Truxal 
1955,  Del  Toro  and  Parker  1960,  Axelby  1960,  Seshadzi  et  al.  1969), 
and  a  standard  transfer  function  is  desired,  this  is  a  design  problem. 
Chen  and  Shieh  (1970)  and  Wakeland  (1976)  have  proposed  analytical 
methods  for  the  compensator  fitting.  However,  their  methods  are 
limited  to  filters  and  compensators  in  which  the  unknown  coefficients 
can  be  solved  by  a  quadratic  equation.  The  method  of  this  paper 
overcomes  this  difficulty. 

(2)  The  transfer  function  obtained  in  this  paper  is  the  function  of  the  original 
system.  When  dominant  frequency-response  data  can  be  obtained 
from  experimental  data  of  a  practical  system  and  the  mathematical 
function  of  the  system  is  desired,  this  is  an  identification  problem. 

(3)  When  the  dominant  frequency-response  data  are  obtained  from  a  given 
high-order  transfer  function  and  various  low -order  approximate  models 
are  required,  this  is  the  model  reduction  problem.  The  reduced  models 
obtained  in  this  paper  have  the  same  selected  dominant  frequency- 
response  data  as  the  original  system.  Thus,  the  design  processes  in  the 
frequency  domain  can  be  greatly  simplified. 


4e2 


2.  Modelling  non-linear  equations 

Given  a  transfer  function  T(s)  of  a  unity  ratio  feedback  closed-loop  system 

T(  ,  60-f ■b1s  +  bg8i+  ...  +  bm8m  _  n{s)  _  G(s) 

a0  +  ats  +  a ^  +  . . .  +  anan  d(s)  1  +  G(s) 

where  n(s)  and  d(s)  are  the  numerator  and  denominator  polynomials,  res¬ 
pectively,  and  at  and  b{  are  constants.  If  the  system  is  a  type  l  system,  the 
open-loop  transfer  function  G(s)  is 

r(  ,  -g(l+c1a  +  caa!i+  +yp)  P(s) 

{)  s‘(l  +dl8+dji*+  ...  +djfl)  q(s)  '  ' 

where  p(s)  and  q(s)  are  the  numerator  and  denominator  polynomials.  K,  l,  c(, 
and  di  are  constants.  AT  is  a  velocity -error  constant  (Kv)  if  1=1. 

The  equations  for  dominant  frequency-response  data  are  : 

(1)  System  type  is  determined  from 

G(jw 0)  =  Re  [<?(>„)]  +  j/m[G!0"ajo)]  at  o>0  =  0  (2  a) 


G(jO)  =  Re[G(jO)] 
T(jO)=~° 

ao 

Re  [<?(iO)]  =  ^t(Cj— dj) 
/J«0’0)]«oo 


for  a  type  0  system 


for  a  type  1  system 


T(jO)  =  —  =  1 
°o 


(2)  Phase  margin  gives 


where 


=  180°  +  L  G(jatc) 


\G(ju>c)\  =  l 

wc  is  the  gain  crossover-frequency. 

(3)  Gain  margin  yields 


where 

L  G[jwr)  =  —  180°  (4  6) 

Wg  is  the  phase  crossover  frequency. 

(4)  Mp  =  |  T(jutv)  [  =  maximum  value  of  the  closed-loop  frequency 


response 

where 


d\T(jw)\ 

du> 


cap  is  the  peak  value  frequency. 
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(5) 

mjcub)i=~ 

(6) 

where  tob  is  the  bandwidth. 

(6) 

IGO'o.^)  )=5-6 

(7  o) 

or 

20  log  |6?(,)'to)(  =  +  15  dB  at  cu  =  cocl 

(7  b) 

and 

lG(Jcoet)  1=0-18 

(7  c) 

or 

20  log  )G(ju>)j  =  —  15  dB  atte  =  o»e2 

(7  d) 

A  set  of  non-linear  equations  can  be  formulated  from  the  basic  definitions 
of  the  assigned  dominant  frequency-response  data  in  (2)-(7).  The  procedures 
can  be  illustrated  by  using  the  following  example.  The  dominant  frequency- 
response  data  in  (2  c),  (3),  and  (4)  are  shown  in  Fig.  1,  which  are  marked  as 
A,  B,  and  C  and  given  as  follows  : 


(1)  Re  [tr(jo>0)]=  — 2-1  and  Im[G(ju>0)] ~  cc  at  o<0  =  0  rad/s 


or 

T(jto 0)  =1  at  £o0  =  0  rad/s 

(8  a) 

(2) 

Re  [G/jaj,)^  —  1*5  at  u>„  =  1-9  rad/s 

(8  6) 

(3) 

LG(ju>n)~  -  180°  at  1-9  rad/s 

(8  c) 

(4) 

^ro  =  180°+  L  G(jw0)  =  5-7°  at  toc  =  3-2  rad/s 

(8  d) 

(5) 

j  G( jwc)  j  =  1  at  o>c  =  3-2  rad/s 

(8e) 

Five  conditions  are  given  in  (8).  Therefore,  various  transfer  functions  with 
five  unknown  coefficients  can  be  constructed.  Assume  that  the  desired 
transfer  function  Td(s)  is 


Ta(s)  = 


bo+b^  +  bys2 
a„  +  ols  +  a^s2  +  ajS3 


(9  °) 


From  the  conditions  in  (8  a),  it  may  be  observed  that  the  system  is  a  type  1 
system.  Therefore  b0=a0.  Also,  to  simplify  the  equation  wo  let  o,=  l. 
Thus,  we  have 


*>)- 


ao+b^+bys* 
a0+a1s  +  aas*  +  s3 


The  corresponding  open-loop  transfer  function  Gd(s)  is 


(9  b) 


«„(»)' 


K(l  -f  CjS  +  CgS*) 
s(l  +dls+dTs1) 


where 


K- 


®o  bt  __  bt  j  at  —  bt 

'  ~  >  ci~~<  c*~  ai~z — ir 

a, -Oj  a0  a0  ot-o, 


(10) 


and  6,= 


®i“&i 
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’  Following  the  basic  definitions  and  the  assigned  data  in  (8)  yields  a  set  of 

non-linear  equations  : 

t 

(1)  The  assignment  in  (8  a),  or  Re  [G(jO)]=  -2-1,  gives 

/i(«».  «i,  bv  bt)=a1bl-b12-a(pt  +  a0bt  +  2l(al-bl)t  =  0  (Ha) 

(2)  The  specification  in  eqn.  (8  b),  or  Re  [G(jiu„)}  ~  -  1*5  at  a>„  =  1-9,  yields 

fi(a0’  ®1)  ®2>  ^lt  ^2)  =  62 )(®0  —  3*6162)  6j (flj  —  bj  —  3*61 ) 

-l*5[3*61(«2-62)1!  +  (a1-61-3*61)1!]  =  0  (11  6) 

(3)  The  condition  in  (8c),  or  LG(ju>„)=  -180°  at  wn=  1-9,  gives 

/s(ao>  ®i»  ®2»  ^i»  ^2) =  3*6161(o1  -  62) 

+  (a0-3*61&2)(al-&1-3*61)  =  0  (lie) 

(4)  The  specification  in  (8  d),  or  <j>m  =  5-l°  at  coc  =  3*2,  yields 

/ 4(^0*  ®i>  ®2»  ^i>  ^2)  lO^h^dj  —  62)  +■  (Oo  —  10*2462)(aj  —  6,  —  10*24) 
-0*319  402  24[(a2-62)(a0-  10*2462) 

10*24)]  =  0  (lid) 

(5)  The  assignment  in  (8  e),  or  )G(jioc))  —1  at  u>c  -3-2,  gives 

fb(ao<  ffln  ®2>  ^i>  ^2)  —  (ao~  10*2462)2+  10*246j2 

—  104*8576(a2  — 62)2  —  10*24(«!  ~6t  —  10*24)2  =  0  (lie) 

Equation  (11)  is  a  set  of  high-order  simultaneous  non-linear  algebraic  equations 
which  are  very  difficult  to  solve.  Considering  the  availability  of  the  computer 
program  package  (IBM  1977)  (called  the  Z  systems)  in  many  digital  computers 
r  for  the  solution  of  non-linear  equations,  the  Newton-Raphson  multidimensional 

method  is  suggested  for  solving  these  equations.  However,  it  is  well  known 
that  the  Newton-Raphson  method  will  only  converge  for  a  small  range  of 
starting  values  or  the  initial  guesses.  A  set  of  good  initial  guesses  must  be 
determined  for  rapid  convergence  of  the  numerical  method.  Four  methods 
are  proposed  for  these  good  initial  guesses. 

3.  The  initial  guess 

It  is  well  known  that  high-order  non-linear  equations  have  many  solutions. 
The  solution  and  the  speed  of  convergence  of  a  numerical  method  depend 
heavily  on  the  initial  guesses  or  the  starting  values.  In  this  paper,  the  Newton- 
Raphson  method  is  suggested  for  solving  the  non-linear  equations.  The 
following  methods,  depending  on  the  applications  of  interest,  are  proposed  for 
good  initial  guesses. 

3.1.  Initial  guess  by  a  synthesis  method 


Suppose  only  the  dominant  frequency-response  data  in  (8)  are  available 
and  an  approximate  transfer  function  Td*(s)  of  the  desired  Ta(s)  in  (9  6)  is 
required.  The  TA*(s)  is 
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where  «,*  and  6,*  are  the  starting  values  of  the  numerical  method.  The  steps 
to  obtain  (12)  are  summarized  as  follows  : 


Step  1.  Determine  a  second-order  approximate  transfer  function  T3*(s) 
using  <£m  =  5-7°  and  u>(.  =  3-2  rad/s  in  (8  d)  and  (8  e).  This  Tt*  (a)  is 


7V(s)  = 


a2  +  2fa»ns  +  cuD2 


(13  a) 


where  f  =  the  damping  ratio  and  ain  =  the  natural  angular  frequency.  Two 
non-linear  equations,  which  are  constructed  from  the  basic  definitions  of  wc 
and  <f>m,  can  be  obtained.  These  non-linear  equations  can  be  converted  into  a 
single  variable  (f  or  ojn)  high-order  equation  from  which  the  roots  can  be 
determined.  Using  this  approach,  we  have  £= 0-0498  and  a>n  =  3-2079.  The 
poles  that  can  be  considered  as  the  dominant  poles  of  a  system  can  be  deter¬ 
mined  from  the  characteristic  equation  in  (13  a).  The  dominant  poles  are 


a, 2=  -f<un  +  —  0-1598  +  ^'3-2039 


Thus,  (13  a)  becomes 


W)  = 


10-2909 

+  0-31 94s +  10-2909 


(13  6) 

(13  c) 


Step  2.  Construct  a  third-order  approximate  transfer  function  T3*(s)  by 
inserting  in  it  a  pole  (a  =  ~p)  and  modifying  the  term  in  the  numerator  of 
T2*(a)  so  that  the  final  value  of  the  Ta*(s)  equals  to  unity,  or 


T3*(s) 


P<*n* 


(a*  +  2  £una  +  wn2)(*  +  P) 


10-2909P 

(*2  +  0-3194s+  10-2909)<a  +  P) 


(13  d) 


The  unknown  constant  P  can  be  easily  determined  by  using  the  condition  in 
(8  6),  or  Re  [f»(jto„))  =  -  1-5  where  «>„  =  1-9.  Thus,  we  have 


P  =  4-5401 


(13  e) 


Step  3.  Establish  another  third-order  approximate  function  T3**{s)  by 
inserting  a  zero  in  (13  d)  with  an  unknown  constant  61*. 


T3**(s) 


6,*  a  +  Pid„! 

(a2  +  2fa>na  +  wn2)(a  +  P) 


6j*  s  +  46-7216 

(a2  +  0-3 194  a  +  10-2909)(a  +  4-5401) 


(13/) 


The  6,*  can  be  determined  by  using  the  condition  in  (2  c)  and  (8  a),  or 
Re  [G(j0)  |  =  -  2- 1 .  The  6,*  is 


Hence,  we  have 


6i*  =  32-4038 

46-7216  +  32-4038* 
46-7216+  1 1-7410* +  4-8595a2  +  a3 


(13  g) 
(13  k) 


Equation  (13  6)  can  be  considered  as  an  approximate  function  of  (12)  by 
a&suming  6g*  =  0.  The  initial  guesses  in  (12)  are  a0*  =  46-7216,  at*  =  11-7410, 
ag*  =  4-8595,  6[*  =  °  .’-4038,  and  6g*  =  0.  Using  these  constants  as  starting 


values  for  the  numerical  method  yields  the  desired  coefficients  in  (9  6),  or 
n0  =  6-378  070,  a,  =  10-462  220,  a2  =1-259  008,  6t  =  20-556  61 ,  and  bt  =  0-243  466. 
The  desired  transfer  function  is 


m  y  6-378  070  +  20-556  61s +  0-243  466s2 
3^“  6-378  070+  10-462  220s  +  1-259  008s2 +  s3  (14) 

The  Newton-Raphson  method  (IBM  1977)  converges  at  the  9th  iteration  with 
the  error  tolerance  of  10_a.  Equation  (14)  has  the  exact  frequency-response 
data  specified  in  (8). 


3.2.  Initial  guess  by  complex-curve  fitting  and  continued  fraction  methods 

The  problem  of  finding  unknown  coefficients  of  a  transfer  function  as  a 
ratio  of  two  frequency-dependent  polynomials  has  been  investigated  by  Levy 
(1959).  His  method  minimizes  the  sum  of  squares  of  the  errors  at  arbitrary 
experimental  points.  We  present  a  simple  method  to  determine  the  approxi¬ 
mate  coefficients  of  a  transfer  function  using  the  real  parts  and  imaginary  parts 
of  available  limited  frequency-response  data.  A  low-order  model  is  often 
determined  because  of  data  limitation.  The  low-order  model  is  then  expanded 
into  a  continued  fraction  of  the  Cauer  second  form  to  obtain  a  set  of  dominant 
quotients.  Then  some  non-dominant  quotients  are  inserted  into  the  continued 
fraction  to  obtain  an  amplified-order  model  (Huang  and  Shieh  1976)  which  is 
the  desired  approximate  transfer  function  for  the  use  of  the  initial  guess. 

Consider  the  transfer  function 


y/j)- *>»  +  *>!«  +  »>«*+  ••• 

1  +o1s+a2s2+  ...  +  a„sn 


(15a) 


where  a(  and  bi  are  unknown  coefficients  to  be  determined.  Substituting 
s  —  j<ok  into  (15  a)  we  have 


r*0A)  = 


(60  -  b2wk2  +  64<V  -  b6wke  +  ...) 

_ +  -  b3wk3  +  -  6t<V  + 

( 1  —  atfnk2  +  atajk*  —  a6oifc*  +  ...) 

+  j  ~  «s  “>**  +  ascuk!‘  ~  ®7  <V  + 


-) 


-) 


=  R(a>lc)+jI(u>k)  =  Rlc+jIk  (15  b) 

when  Rk  and  lk  are  the  given  real  and  imaginary  parts  of  the  T*(s)  at  the 
available  frequencies  u>k.  Multiplying  both  sides  of  (15  6)  by  the  common 
denominator  and  separating  the  real  and  imaginary  parts,  and  also  equating  the 
respective  real  and  imaginary  parts,  yields 

60  -  b^Mk%  +  bAu)k*  -  6«<u(.®  +  ...  +alIkwk  +  a2Rkwk2 

~  a3^kwk3  ~  a^Rkwk*  d"  ■•■~Rk  (15  c) 

and 

b\o>k  —  63id*3  +  bku>ks - bjw*7  +  ...  -atRkwk  +  atIk wk2 

+  aiRkwk3-aiIkwk*+  ... -Ik  (15  d) 
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In  matrix  form,  (15  c)  becomes 


1 

“"l* 

wi* 

—  tO,6 

■ff,to,2 

-/,a>i3 

—  i?,to,4  . 

p.  i 

r^ii 

1 

—  to,2 

to24 

—  aq8 

-/gtOj3 

—  ifjCUg4  . 

** 

1 

-a>32 

"3* 

-a»36  . 

i?3a)32 

-h^a3 

-R3cj3*  . 

^4 

ai 

(15  e) 

«2 

■ 

1 

-a>3 

“z4 

~o>z®  ' 

4^z 

Rxu>x2 

"h 

3 

• 

1 

-Rxwx*  . 

.Rx_ 

where  x  =  n  +  to/ 2  +  1  if  to  is  even  and  x  —  n  +  (to  +  1  )/2  if  to  is  odd. 

Substituting  a,-  obtained  in  (15  e)  into  (15  d),  we  have  another  matrix 
equation  to  solve  for  bit  i  =  1 ,  3,  5 . 


'to, 

-to,3 

to,3 

-to,7 

“ 

r6‘  i 

to2 

-to23 

to26 

-to27 

*3 

to3 

to33 

-to37 

to„3 

-to„7 

.5,. 

((«„/ ,to,°  +  alR1w11)  -  (a2J,a>,2  +  a3R1u>13)  +  ...) 
2OJ2°  a^R2u)^ )  (o2/2to22  +  a3R2w23)  +  . . . ) 


L««oW + ai^voiv1 )  —  (a2I  ua>v2 + a3Rvcov3)  +  ...)J 


(15/) 


where  tofc°=l,  o0=l  ;  fc  =  w  and  y  =  (m+l)/2  if  m  =  odd  ;  &  =  to— 1  and  y  = 
to/2  if  to  =  even.  In  this  example,  the  available  data  are 


Since  only  three  values  are  available,  the  approximate  function  Tt*(s)  is 


f,*(*)- 


bp  +  bjs 
1  ■(■  fljfl  “I- 


(17o) 


Substituting  the  data  at  to„  and  to2,  and  to3  in  (16)  into  (15  e)  yields  fc0  =  1, 
a,  =  0-0388,  and  oa  =  0-1839.  Then  substituting  o(  and  the  data  at  to,  into 


(15  /)  gives  6t  =  2-8907.  Because  the  desired  approximate  function  in  (12)  is  a 
third-order  function,  T2*(s)  should  be  amplified  by  using  the  continued  fraction 
method  (Huang  and  Shieh  1976)  as  follows. 

T2*(s)  is  first  expanded  into  a  continued  fraction  of  the  t'auer  second  form 
to  obtain  a  set  of  dominant  quotients  :  hl  =  \,  h2=  —0-3507,  h3  =  —0-9651,  and 
A4  =  16-0725.  Then  the  order  of  T2*{s)  is  amplified  to  the  third  order  by 
inserting  non -dominant  quotients  A5=  100  and  A6  =  0-1,  or 


1  +  2-8907s 


1 


+  0-0388s  + 0-1 839s2  ,  s 

A,  4 - 


h j  +  - 


he  4  - 


^3  +  T- 


he  4-  - 


K 


he+- 


4  - 


i  s 

hi  4  — 


=  T.*(. 


54-3885  4  162-6914s  4  15-8219*2 
'  54-3885  4  7-5839, s  4  1 0-2 146s2  4  .s3 


h» 

(17  b) 


Huang  and  Shieh  (1976)  have  shown  that  the  amplified-order  model  is  a 
good  approximation  of  the  original  low-order  model  if  the  inserted  positive 
quotients  h( p  1  and  hi+t4  1  where  i  is  an  odd  number.  Using  the  coefficients 
in  (17  b)  as  initial  guesses  we  have  the  desired  coefficients  in  (14)  at  the  15th 
iteration  (IBM  1977)  with  the  error  tolerance  of  10  ® 

If  much  experimental  frequency -response  data,  including  the  dominant 
data  of  a  system,  is  available  and  the  transfer  function  of  the  original  system  is 
required,  this  is  an  identification  problem.  In  this  case,  a  set  of  non-linear 
equations,  based  on  the  basic  definitions  of  the  dominant  data,  can  be  con¬ 
structed  and  can  be  solved  by  the  Ncwton-Raphson  method.  The  initial  guess 
can  be  determined  by  using  the  dominant  data  and  others  in  (15).  Since  main- 
data  are  available,  a  high-order  approximate  transfer  function  can  be  deter¬ 
mined.  Therefore,  the  use  of  the  continued  fraction  method  (Huang  and 
Shieh  1976)  is  not  necessary. 

When  a  high-order  transfer  function  of  a  system  is  given  and  various 
reduced-order  transfer  functions  are  required,  this  is  a  model  reduction  problem. 
In  the  frequency  domain ,  numerous  methods  (Chen  and  Shieh  1969,  .Shieh 
and  Goldman  1974,  Hutton  and  Friedland  1975,  Sharnash  1975,  Lai  and  Van 
Valkenburg  1976)  have  been  proposed  for  model  reduction.  The  continued 
fraction  methods  (Chen  and  Shieh  1969,  Shieh  and  Goldman  1974),  the  Routh 
approximation  method  (Hutton  and  Friedland  1975),  the  time-moment 
matching  method  (Sharnash  1975),  and  the  frequency-moment  matching 
method  (Lai  and  Van  Valkenburg  1976)  are  the  typical  examples.  These 
methods  have  been  critically  compared  by  Decoster  and  Cauwenberghe  (1976). 
The  new  method  presented  in  this  paper  can  be  used  to  obtain  the  reduced- 
order  models  which  have  the  exact  dominant  frequency-response  data  as  those 
of  the  original  one.  This  method  can  lie  called  a  dominant  frequency-response 
data  matching  method.  The  procedure  is  as  follows. 
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Step  1.  Plot  the  frequency-response  eurves  to  determine  the  data  at  the 
dominant  frequencies  u>0,  io„,  to,.,  to,.,,  to,.2,  to,,,  and  to,,. 

Step  2.  Formulate  a  low-order  model  with  unknown  coefficients,  and 
write  a  set  of  non-linear  equations  based  on  the  basic  definitions  of  the  data 
at  dominant  frequencies. 

Step  3.  Determine  a  set  of  good  starting  values  by  using  the  synthesis 
method  or  the  complex  curve  fitting  method,  and  solve  the  non-linear  equation 
by  using  the  Newton-Raphson  method.  Thus,  reduced-order  models  can  be 
determined.  Comparing  the  reduced-order  models  obtained  from  the  proposed 
method  with  those  of  the  existing  methods  (Chen  and  Shieli  1909,  Shieh  and 
Goldman  1974,  Hutton  and  Friedman  1975,  Shamash  1975,  Lai  and  Van 
Valkenburg  1970),  we  observe  that  the  model  obtained  in  this  paper  is  superior 
to  existing  methods  in  that  the  reduced  model  has  the  exact  dominant  frequency 
response  as  the  original.  As  a  result,  an  engineer  can  design  a  control  system 
more  efficiently  in  the  frequency  domain. 

Since  the  original  high-order  transfer  function  is  available,  an  existing 
method  (Chen  and  Shieh  1969)  can  be  applied  and  modified  to  obtain  an 
approximate  transfer  function  for  the  determination  of  the  initial  guess.  Two 
additional  methods  for  initial  guess  determination  are  as  follows. 

(3)  Initial  guess  by  a  continued  fraction  method  (Chen  and  Shieh  1969). 

Consider  the  high-order  transfer  function  in  (1  a).  The  function  can  be 
expanded  into  a  continued  fraction  and  various  reduced  models  obtained  by 
discarding  some  of  the  quotients,  or 


j .  =  6,  +  y+  +bmsm  =  n(s) 

«0  +  al«+  •••  +«n5"  d(9) 


(18  a) 


1 


hy  + 


1 


A, 

hyhty  +  9 


I 


A,+ 


*.+ 


,  « 
A»+r 


h 2^1  ■yhy  -f-  (// g  4"  hyjs 

^  (hyh2  4"  AjA,  4“  h^hy^s  4"  s 


(18  b) 

(18c) 


(18  d) 


Using  the  coefficients  of  the  approximate  model  in  (18)  as  the  initial  guess  for 
the  numerical  method,  we  have  the  desired  reduced  model.  However,  the 
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1 


» 


approximate  model  in  (18)  may  be  unstable  even  if  the  original  system  is  stable. 
The  continued  fraction  method  (Chen  and  Shieh  1969)  can  be  modified  by  the 
following  new  method. 

(4)  Initial  guess  by  a  mixed  method  of  the  continued  fraction  approach  and 
Gustafson’s  (1965)  method. 

Assume  the  reduced  model  of  the  original  system  in  (18  a)  is 


T  *h)  bo*  +  bi*8+  +6p-i*'sP~1  n*(8) 

v  a0*  +  a1*s+  ...  +ap*s1>  d*(s)’ 


ap  =  l 


(19  a) 


A  matrix  equation  (Chen  and  Shieh  1970)  can  be  constructed  from  the  dominant 
quotients  h{,  i  =  1,  2,  ....  p,  obtained  in  (18  b)  and  the  unknown  coefficients  at* 
and  bf*  in  (19  a)  as 


where 


[b]  =  [ff][a] 

(19  6) 

[a]T  =  [a0*,  a,*,  ... ,  ap_, *] 

(19c) 

[bf  =  [b0*  b,*,...,  &„_,*] 

(19  d) 

(19  e) 

where  T  designates  transpose, 


[H2]  = 


*  K 

0 

0  . 

0 

0  1 

r  i 

0 

0  . 

0 

0  1 

1 

^2  0 

0 

0 

0 

K 

0  . 

0 

0 

0 

l  h3 

0 

0 

0 

i 

^2 

0 

0 

0 

0 

1 

0 

0 

0 

0 

j 

0 

0 

0 

0 

0  . 

1 

K. 

0 

0 

0  . 

1 

hp- 1  . 

[#.]  = 


1 

0 

0  . 

0 

0  n 

r  i 

0 

0  . 

0 

0  1 

0 

K 

0  . 

0 

0 

0 

1 

0  . 

0 

0 

0 

1 

h3  . 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0  . 

1 

V 

0 

0 

0  . 

1 

^p-1  . 

1  0  0 


0  0 


0 

1 

0 

.  0 

0 

0 

0 

1 

.  0 

0 

0 

0 

0 

.  0 

0 

0 

0 

0 

.  0 

*1 . 

■  1 

0 

0 

.  0 

0  ' 

0 

1 

0 

0 

0 

0 

0 

1 

.  0 

0 

0 

0 

0 

.  0 

0 

0 

0 

0 

.  0 

A,. 

The  a{*  in  (19  c)  can  be  determined  from  the  coefficients  of  the  polynomial  that 
is  obtained  from  the  product  of  the  dominant  eigenvalues  of  the  d(s)  in  (18  a). 
When  the  dominant  poles  of  d(s)  cannot  be  clearly  identified  or  the  poles  of 
d(a)  are  not  available,  the  paper  and  pencil  method  suggested  by  Gustafson 
(1965)  can  be  applied  to  construct  the  d*(s)  or  to  determine  a(*  in  (19  c).  Then, 
substituting  the  at*  into  (19  b)  yields  the  required  n*(s)  or  b(*  in  (19  a).  The 
steps  determine  the  d*(s)  are  shown  as  follows. 
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Step  1.  Construct  a  Routh  (1877)  array  using  the  coefficients  ai  of  d(s) 
and  the  Routh  algorithm.  The  a(  are  expressed  by  double-subscripted  notation 
atJ  for  obtaining  the  general  algorithm.  The  Routh  array  is 


In  general  — ®<_2,/+i— 'y<-2®<— i,j+i  ;  i—  1,  2,  ...,  j  —  3,  4,  ... 

Yi=ai,ilai+i,i  (20fc) 

Step  2.  Construct  various  approximate  low-order  polynomials  dt*(s)  from 
the  last  row  and  the  next  to  last  row,  and  so  on  in  the  Routh  array. 

For  example,  the  ith  order  approximate  equations  are 

^x*(s)  =  «»,xs  +  «n+i.i=0B,i5  +  ®o  =  0  whent=l  (20c) 

d**(a)  =  a„_i,  j«*  +  ®B>  4-  «B_i,  2  =  an-i,  i®2  +  an,  +  a„  =  0  whent  =  2  (20  tf) 

and 

=  ®n-2,  Is*  +  ®b-1, 11*2  +  an-2,  2s  +  ®b-1,  2 

=  «B-s, x«*  +  «»-i,  X®2  +  °n-2, 2 s  +  ®0  =  0  when  *  =  3  (20  e) 

Since  the  original  system  is  asymptotically  stable,  all  are  positive  values. 
The  approximate  polynomials  dt*(s)  are  always  the  Hurwitz  polynomials. 
Moreover,  Gustafson  (1965)  has  shown  that  relationships  exist  between  the 
coefficients  of  di*(s)  and  the  time-domain  moments.  The  normalized  poly¬ 
nomials  can  be  determined  by  dividing  each  coefficient  in  dt*(s)  by  the  coeffi¬ 
cient  of  the  highest  order  term  in  s.  The  approximate  transfer  function  Tp*(s) 


•i 
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in  (19  a)  can  be  considered  as  a  reduced-order  model  of  the  original  high-order 
1  system.  In  this  paper,  we  use  it  as  the  initial  guess  for  the  numerical  method 

for  determining  the  reduced  ofder  mode)  that  has  the  exact  dominant  frequency- 
response  data  as  the  original  system. 


4.  An  illustrative  example 

Consider  the  unit  ratio  feedback  closed-loop  transfer  function  of  a  stabilized 
real  missile  system  (Bosley  1977) 


where 


and 


T(s)  =  '  ”  1 

■  ■  o  / 

a0  +  axs  + 

...  +ansn 

a0  =  8-802  158  509  x  1018, 

a,  =  2-41 9  047  424  x  1019 

a2  =  2-91 1  920  56  x  1018, 

a3  =  2-420  405  431  x  1018 

<z4  =  6-667  397  031  x  1018, 

a5  =  9-749  923  212  x  1014 

ag  =  9-360  329  977  x  1012, 

a,  =  6-231  675  318  x  1010 

ag=  2-976  950  696  x  108, 

a9  =  9-316  239  04  x  10s 

a10=  1-923  554  x  108, 

au=  1 

k'  =  1-494  523  312  x  1011 

6'0  =  5-889  609  375  x  107, 

b\  =  3-084  598  703  x  108 

f>'2=  1-958  045  299  xlO7, 

6'3  =  3-357  065  095  x  10s 

6'4=  1-715  193  3  x  10s, 

II 

order  and  the  third  order  reduced-order  models  which 

(21  a) 


of  the  dominant  frequency-response  data  of  the  original  system  are  required. 
The  open-loop  transfer  function  G(s)  of  the  system  is 


G(s)  = 


k(e0  +  els+  ...  +e5s5) 
s(?o  +  Sh»+  +g10s10) 


(216) 


where 


and 


5T0  =  -2-190  952  724  6  x  1019, 

gt=  -  1-442  378  55xl018 

gt  =  2-370  233  311  x  lO18, 

gr3  =  6-641  763  067  x  1016 

g4  =  9-748  428  689  x  1014, 

g&=  9-360  329  977  x  1012 

gt  =  6-231  675  318  x  1010, 

g7=  2-976  950  696  x  108 

gs  =  9-316  239  04xl05, 

gr9=  1-923  554x  10* 

9 10  —  1 

k=  1-494  523  312  x  10u 

e0  =  5-889  609  375  x  107, 

e,  =  3-084  598  703  x  108 

e2=  1-958  045  299  x  107, 

e3  =  3-357  065  095  x  10* 

e4=  1-715  193  3 x  10*, 

*5=1 
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Note  that  (>'{s)  is  u  non-minimum  phase  function  ;  its  Nyquist  plot  is  shown  in 
Fig.  1.  The  dominant  frequenc 'y -response  data  are  chosen  and  given  in  (8). 
The  set  of  non-linear  equations  are  shown  in  (11).  The  initial  guesses  show:1 
in  (13  /i)  and  (17  h)  yields  the  required  third-order  reduced  model  in  (14),  or 


0-243  466.-r  +  20-551;  <11* +  0-378  <>7 
*3  +  1  -250  008.v2  +  1 0-402  22s  +  (>-378  07 


(22  a) 


If  the  continued  fraction  method  (Chen  and  Shieh  1000)  in  (18)  is  used,  the 
approximate  reduced  model  is 

0-0020*-  +  1 0-4092*-  +  3-7370 

T.  ,*(#)  =  - - (22  6) 

31  .v3  +  0-0488, >'3  +  10-1 001* +  3-7370  '  ' 


Using  the  coefficients  in  (22  b)  as  starting  values  for  solving  the  non-linear 
equations  in  (11)  yields  the  desired  coefficients  in  (22  a)  at  the  eighth  iteration 
(IB3I  1077)  with  the  error  tolerance  of  10'  .  If  the  mixed  method  in  (10) 
and  (20)  is  used,  the  normalized  approximate  denominator  in  (20  e)  is 

,/,*(#)=«» +  0-9524#*  +  10- 1924#  +  3-7455  (22  c ) 

The  ii3*(s)  obtained  from  (19)  is 

n3*(s)  =  0-7000*--  +  19-5155#  +  3-7455  (22  </) 

The  approximate  transfer  function  by  the  mixed  method  is 

n,  tki  4  0-7000#*+ 19-5155#  +  3-7455  /M  ^ 

3m  ^  =  #3  +  0-9524#*  +  1 0- 1 924#  +  3-7455  ’ 


If  the  coefficients  in  (22  e)  are  used  as  starting  values,  the  Newton-Haphson 
method  ( I fi.il  1077)  will  converge  to  the  desired  solution  in  (22  n)  at  the  eighth 
iteration  with  the  error  tolerance  of  10  s.  The  unit  stop  response  curves  of 
various  reduced  models  and  the  original  system  are  compared  in  Fig.  2.  All 
three  reduced-order  models  give  very  satisfactory  approximate  time  response 
curves.  However,  only  the  7\,*(.v)  in  (22  a),  which  uses  the  method  of  dominant 
frequency-response  data  matching,  has  the  exact  dominant  frequency -response 
data  as  the  original  system. 

If  oj,.  =  3-2  rad/s,  <f>m  =  5-7J  and  Re  K>'(/0)|=  —2-1  are  chosen  as  the  domi¬ 
nant  data,  the  second-order  reduced  model  obtained  by  the  proposed  method  is 


3-339  51 7#  +  9-224  24 
s2  +  0-302  805s  +  9-224  24 


(23  a) 


The  approximate  reduced  models  by  the  continued  fraction  method  and  the 
mixed  method  are  : 


m  24*7081.?  -f  4*8122 

*  '  #*  +  12-82(1 1«  +  4-8 122 

(23  b) 

m  15-3818#  +  3-9328 

*"■  #*  +  6-5726#  +  3-9328 

(23  c) 
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Figure  2.  Time  responses  or  original  and  third-order  reduced  models. 


The  unit-step  time  response  curves  of  various  reduced-order  models  T3*(s), 
T2*{s),  T2c*(s),  and  T2m*{s)  are  compared  in  Fig.  3.  It  is  observed  that  7’2*(.s) 
gives  better  approximation  in  the  transient  response  than  Ttl.*(s)  and  Tin*(s). 


Figure  3.  Time  responses  of  third-  and  second-order  reduced  models. 
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5.  Conclusion 

A  basic  method  has  been  developed  for  modelling  transfer  function  using 
dominant  frequency-response  data.  When  the  specifications  of  the  design 
goals  of  a  control  system  are  assigned,  the  proposed  method  gives  the  standard 
transfer  function.  Tims,  the  design  processes  in  the  frequency  domain  can  be 
significantly  simplified.  When  the  experimental  frequency-response  data  of  a 
system  are  available,  the  proposed  method  can  be  used  to  identify  the  transfer 
function  of  the  original  system.  Also,  if  a  high-order  transfer  function  is 
given,  various  low-order  models  can  be  determined.  The  reduced  models  have 
the  same  dominant  characteristics  of  the  original  system.  Four  methods  have 
been  proposed  for  estimating  the  good  starting  values  for  the  solution  of  non¬ 
linear  equations.  The  new  dominant  frequency-response  data  matching 
method,  and  the  new  mixed  method  that  has  the  advantages  of  both  continued 
fraction  method  of  Chen  and  Shieli  (1909)  and  the  paper  and  pencil  method  of 
Gustafson  (1905)  have  been  developed  for  model  reduction. 
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